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Numerical  Calculation  of  Steady  Two-dimensional 
Axi-symmetrical  Turbulent  Boundary  Layer 
of  a  Compressible  Fluid 


Wang  Ying-shih 


(Institute  of  Mechanics,  Academia  Sinica) 


This  paper  employed  a  momentum  equation 
expressed  by  the  stress  tensor  to  derive  the  momen¬ 
tum  equation  of  a  non-rotational  steady  axi- 
symmetrical  turbulent  boundary  layer  of  a  compressible 
fluid  in  the  axi-symmetrical  coordinate  system.  It 
was  found  that  If  the  flow  in  the  duct  was  to  be 
treated  as  viscous  fluid  everywhere  then  the  curva¬ 
ture  of  the  duct  and  the  rate  of  variation  of  the 
curvature  along  the  direction  of  the  flow  cannot  be 
too  large.  This  has  not  been  proven  in  any  existing 
literature . 


Based  on  the  set  of  equations  proposed  in  this 
paper,  we  compiled  a  program  and  computed  an  example. 
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In  the  actual  fluid  flow  process,  especially  under  heat 
transfer  and  mass  transfer  conditions,  the  effect  of  viscosity 
is  very  important.  Furthermore,  most  of  the  viscous  flow 
process  is  a  turbulent  flow.  Therefore,  in  order  to  solve 
this  type  of  a  flow  process  perfectly  it  must  begin  with  the 
solution  of  the  turbulent  viscous  flow  equation  set.  However, 
the  work  load  involved  in  obtaining  solution  to  this  set  of 
equations  is  rather  large.  Until  the  late  sixties  with  the 
development  of  computational  methods  and  computers  it  was  then 
possible  to  work  in  this  area. 

The  viscous  flow  problem  can  be  divided  into  viscous  flow 

with  return  flow  and  the  boundary  layer  flow  without  return 

flow.  These  two  types  of  flow  have  different  types  of  basic 

equations.  Their  solution  and  the  treatment  of  boundary  condi- 

T  1-^1 

tions  are  not  the  same1-  JJ.  This  paper  mainly  discussed  the 
calculation  of  fluid  field  of  the  turbulent  boundary  layer  of  a 
compressible  fluid  in  a  two-dimensional  axi-symmetrical  coordi¬ 
nate  system. 

I.  Basic  Equations  and  Discussion  of  the  Content 

When  a  compressible  fluid  undergoes  non-rotational  axi- 
symmetrical  flow  in  a  axi-symmetrical  coordinate  system  as 
shown  in  Figure  1,  the  boundary  layer  kinetic  equations  under 
steady  state  are  shown  as  follows  (detailed  derivation  see 
Appendix  1): 


Figure  1.  Axi-symmetrical  coordinate  system.  Key:  1. 
boundary,  2.  axis  of  symmetry. 
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mass  transfer  equation 
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energy  equation  W  +  .  _i.  JL  iljk  +  V  w,y,  -  nr)  r}  (4) 
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then 


If  the  fluid  satisfies  the  conditions  of  a  gas  completely. 


—  —  RT  - 
P 


If  we  neglect  the  variation  of  C  for  each  component  and 
2  P 

also  neglect  v  ,  then  the  stationary  enthalpy  can  be  defined 


A*  *»  C,T  +  ^  (6) 

l 

The  above  equations  are  applicably  in  the  laminar  flow 
region  of  the  boundary  layer  as  well  as  in  the  turbulent  region 
of  the  boundary  layer.  When  it  is  used  for  turbulent  boundary 
layer,  all  the  physical  quantities  should  be  expressed  using 
the  time  average  values. 


At  the  moment  let  us  assume  that  the  momentum,  mass  trans¬ 
fer,  and  heat  transfer  respectively  obey  Newton's  law,  Fisk's 
law,  and  Fourier's  law  in  the  turbulent  boundary  layer  similar 
to  what  happens  in  the  laminar  flow  boundary  layer.  The 
exchange  coefficient  in  the  entire  boundary  layer  can  be 
expressed  using  the  effective  values,  i.e. 
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(7) 


(8) 
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\  <JkM  ' 

\dy  > 

In  the  layer  flow  region  <t*m,  and  <*»•«  correspond 

to  the  viscosity,  Prandtl  number,  and  Schmidt  number  of  the 
fluid,  respectively.  The  effect  viscosity  inside  the  turbulent 
flow  region  is  defined  by  adopting  the  mixing  length  concept 
of  Prandtl  which  is: 

”  pi 1  |  ~  |  (10) 

and  ahM  and  <*iM  are  the  effective  turbulent  flow  Prandtl 
number  and  the  effect  turbulent  flow  Schmidt  number,  respect¬ 
ively. 

Equations  (2),  (3),  and  (4)  are  of  the  conservation  equa¬ 
tion  type.  The  left  hand  side  of  the  equal  sign  is  the  convec¬ 
tion  term.  The  first  term  on  the  right  is  the  diffusion  term 
and  the  second  term  is  the  source.  If  the  source  term  is  known 
or  can  be  obtained  through  another  expression,  then  the  above 
ten  equations  can  give  solution  to  the  ten  unknowns  p,  T, h*,  u,  v, 
•*>,*,  Ji,  ]*>  t*u.  .  During  the  solution  seeking  process  the 
values  for  I,  <riM,  <rt.u  can  be  obtained  from  certain  expres- 

r  p  n 

sions  or  experimental  dataL  . 

In  order  to  obtain  the  solution  more  conveniently,  let 
us  transform  the  above  equations  into  the  von  Mises  coordinate 
system  coordinate  system)  and  bring  in  the  concept  of  the 

flow  function. 


d<]>  dd) 

af  "  pru‘  d7  f 


Because 


(£).-(£).—  (£). 

{%).  -"'“Hi). 

Equations  (2),  (3),  and  (4)  then  become 


S-“  -  JL  (r,) 

dx  dip  pu  dx 

(2*  ) 

5a.  -  a,)  +  & 

dx  dip 

(3’) 

dh * 

dx 

(4  ’  ) 

Because  equations  (7),  (8),  and  (9)  can  be  transformed 
into  partial  derivatives  with  respect  to  t,  J  ^ ,  and  Jh 
respectively,  therefore  equations  (2’),  (3*),  and  (4’)  become 
the  classical  parabolic  curve  equation  with  a  source.  Hence¬ 
forth,  it  is  in  principle  possible  to  adopt  numerical  methods 
which  are  suitable  for  parabolic  equations  to  process  the  above 
set  of  equations. 

Before  proceeding  further  in  seeking  for  the  solutions  of 
the  above  equations,  we  are  going  to  carry  out  a  discussion  on 
the  above  mentioned  momentum  equation. 


(12) 

(13) 


In  the  derivation  of  momentum  equation  in  the  axi- 
symmetrical  coordinate  system  (see  Appendix  I),  the  first 
obtained  equation  Is 


du  .  du 


dr  +  co*  9  r  _  dp 
dy  *  dx 


(1-4) 


Only  general  boundary  layer  assumptions  are  made  in  the 
deviation  process.  The  applicable  region  of  equation  (1-4)  is 
relatively  wider.  It  can  be  applied  not  only  to  the  boundary 
^ver  .ow  of  a  revolving  surface  on  the  basis  of  the  axi- 
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symmetrical  coordinates,  but  also  to  the  jet  boundary  layer 
flow  on  the  basis  of  cylindrical  coordinates. 


As  for  the  analysis  of  axi-symmetrical  jet  flow  it  is 
generally  better  to  use  the  cylindrical  coordinate  system.  At 
this  time  the  thickness  of  the  boundary  layer  and  the  radial 
length  belong  to  the  same  order  of  magnitude.  Therefore,  the 
momentum  equation  of  the  boundary  layer  of  a  steady  two- 
dimensional  axi-symmetrical  jet  flow  can  be  derived 


du  ,  du  1  3  t  \  dp 
-  +  pP  — ■»  -  -  CrO  -  -f- 


where  u  and  v  are  the  velocity  components  along  the  x  and  r 
directions,  respectively. 


Comparing  equation  (1-4)  with  equation  (14)  and  also  referr¬ 
ing  to  Figure  1,  we  can  see  that  when  »-*o  and  r,  —  o  the 
orthogonal  axi-symmetrical  curve  coordinate  system  is  reduced 
to  a  cylindrical  coordinate  system.  Equation  (1-4)  becomes 
equation  (14).  Earlier  Mangier,  in  his  treatment  of  the  laminar 
boundary  layer  problems,  did  not  keep  the  second  term  on  the 
right  hand  side  of  equation  (1-4)  through  order  of  magnitude 
comparison.  Mathematically,  it  transformed  the  boundary  layer 
problem  of  a  revolved  body  into  a  planar  boundary  layer  problem 
for  processing  in  order  to  simplify  the  procedure  necessary  to 
obtain  the  solution  .  Presently,  because  of  the  use  of  numeri¬ 
cal  computations  and  through  the  establishment  of  computational 
programming  based  on  the  equations  stated  above,  we  not  only 
can  solve  the  boundary  layer  problem  of  a  revolved  body  but  also 
can  deal  with  the  boundary  layer  problem  of  the  axi-symmetrical 
jet  flow  and  the  boundary  layer  problem  of  the  planar  flow. 


In  order  to  raise  the  accuracy  of  calculation  in  the  various 
regions  in  the  boundary  layer,  the  choice  of  the  difference  net- 


work  before  the  establishment  of  the  difference  equations  is 
very  important.  If  we  use  the  simple  *-<!»  coordinate  system  to 
establish  the  network  then  in  the  region  where  x  is  small  it  is 
not  possible  to  clearly  express  the  variation  of  the  flow 
field  due  to  the  small  number  of  network  nodal  points.  If  for 
this  reason  the  coordinate  lines  for  <l>  value  becomes  much  more 
dense  in  order  to  have  more  nodal  points  in  the  x  direction  in 
the  area  where  x  is  small,  then  there  are  large  number  of 
unnecessary  nodal  points  in  the  area  where  x  is  large. 

In  overcoming  the  above  difficulty,  some  authors  used  the 
x  —  h/u„  or  x  —  y/ys  coordinate  system.  In  this  paper  we  used 
the  Patankar-Spalding  coordinate  system  mentioned  in  Reference 
[2]  which  is  the  x-oj  coordinate  system.  This  type  of  method 
which  transforms  the  ordinate  of  the  orthogonal  coordinate  sys¬ 
tem  into  a  relative  value  is  to  arrange  so  that  the  number  of 
network  nodal  point  to  be  the  same  during  every  step  as  the 
computation  progresses  in  order  to  improve  the  accuracy  of  the 
computation  in  the  area  where  x  is  small.  Based  cn  Reference  [2] 

u  „  (15) 

<l>u  — 


Figure  2.  Diagram  for  Coordinate  Transformation. 

Key:  1.  Constant. 

Here  the  coordinate  relationships  between  x.l‘,  and  i*>  are 
shown  in  Figure  2.  From  this  we  get  the  following  equations 


(d_\  _ 

f_a_\ 

—  u>  —  (4>  t  —  <£/)] 

(_d_\  + 

dx 

dx 

\dx  J* 

\5x  /■ 

<  Qia  J  * 

<Ps~  4>i  J 
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where 


d±i 

dx 


d-&  -  -n*’,'. 


•  —rtm  i 


fi  ft 

ihj  and  ihg  represent  the  rate  of  mass  transfer  across  the  inner 
and  outer  boundary,  respectively. 


Using  equations  (16)  and  (17)  we  can  transform  the  set  of 
conservation  equations  (2'),  (3')»  and  (4')  into  the  following 
general  form: 


where 


90 

dx 


+  (a  +  bto)  — 
0(0 


a  ■«  r,ih',' /(.<!>,  —  4>i) 
b  —  (r*rf»*  —  r,m',')/C<bs  —  <Pi) 

e  —  r*puf ~ 


(18) 


The  source  term,  however,  has  different  content  which 
varies  with  the  characteristic  of  the  conservation  equation1. 


Characteristic  Content 

of  the  of 

Equation  Variables  term  d 

l  dp 

momentum  equation  u  — — 

mass  transfer  m  JL_ 

equation  j  G>«) 


energy  equation  h* 


9(«,/2)t 
d(o  J 


1The  source  term  of  the  energy  equation  is  derived  after 
letting  <*/.n - <7,.„ 


Equation  18  is  still  a  parabolic  function  and  its  numeri¬ 
cal  solution  is  obtained  through  the  establishment  of  a  differ¬ 
ence  equation  described  in  the  following  section. 

II.  The  Establishment  and  Solution  of 
the  Difference  Equation 

Let  us  take  a  finite  controlled  volume  in  the  flow  field 
as  shown  in  Figure  3  and  then  carry  out  integration  within 
the  limits  of  the  defined  finite  control  volume  to  establish 
the  difference  equation.  The  use  of  this  method  to  establish 
the  difference  equation  has  apparent  advantages  over  the  use  of 
the  Taylor  series  expansion  method  whether  from  the  point  of 
view  of  satisfying  the  physical  concepts  of  the  conservation 
equations  or  for  the  prevention  of  serious  mathematical  error^^' 

Based  on  Appendix  II,  the  difference  equation  is  as 
follows: 


£i“Ao+  +  gl®U  +  g)Qo~  +  g*  “  £»( &D+  ~~  Qd)  ~~  gt (<Po  —  o~ ) 

+'-+ (*).<♦“-♦■> 

After  simplification 


(19) 


<Pu  ™  +  C 


(20) 


where 


A - iLZJb - 

gt  +  gi  +  gt  —  (.dA/W), 

B _ lL=Jl - 

gt  +  £»  +  £«  —  (drf/9$), 

c  -  d‘  -  -  ii 

gt  A-  gi  A-  gt  —  ( dd/d <P). 


For  the  conservation  of  momentum  equation,  due  to  the  differ- 
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ent  contents  from  treatment  of  the  source  term  (see  Appendix 
II),  the  difference  equation  corresponding  to  equation  (20) 
is: 

uu  —  Atr0+  +  B„uu-  +  C,  (21 ) 

where  v  r-<.  +  s, 

“  gi  +  +  gt  —  s> 

B  g>-i>  +  S> 

‘  *»  +  *+*-  S, 

C. - - 

gi  +  gi  +  gi,  —  Si 


Based  on  the  derivation  in  the  Appendix  it  was  found  the 
equation  (19)  is  the  six  point  unexplicit  difference  form  which 
has  the  characteristic  of  stability  for  any  progressing  steps. 


The  coefficients  of  the  linear  algebraic  equations  repre¬ 
sented  by  equation  (20)  form  a  three  diagonal  linear  matrix  and 
solution  can  be  found  using  an  iterative  method.  For  that 
equation  (20)  is  simplified  further  as: 

<P,  -  A'Am  +  b;  (22) 


where 


A\ 


1  -  B,AUX  * 


b; 


BiB!,,  ±  C, 

l-  BiA'i-x  * 


A\ 


A„  Bi 


BA  +  c, 


Because  and  are  determined  by  the  boundary  condi¬ 
tions,  when  the  <J>  1  and  values  are  not  given  on  the 

boundary,  then  their  coefficients  are  zero.  Here  we  theoreti¬ 
cally  solved  the  problem  in  seeking  for  a  solution. 
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Figure  3.  Controlled  Volume 

When  the  nodal  points  are  chosen  to  be  more  dense,  then 
the  linear  distribution  of  $  value  assumption  in  Figure  3  (On 
the  4>-  <i>  plane)  has  already  met  the  accuracy  requirement  for 
all  the  points  in  the  flow  field.  But  because  the  variation 
of  4>  value  is  steeper  near  the  boundary,  if  we  still  use  the 
linear  variation  of  the  true  value  on  the  boundary,  distortion 
will  occur.  Especially  when  the  values  of  relative  flux  and 
transfer  rate  of  a  physical  quantity  on  the  boundary  are  to  be 
determined,  it  would  cause  a  large  error.  Because  when  these 
values  are  determined,  it  is  necessary  to  use  the  gradient  of  <j> 
on  the  boundary  and  this  gradient  can  not  be  replaced  by  the 
simply  linear  variation  between  the  true  $  value  on  the 
boundary  and  the  <5^  value  at  its  neighboring  nodal  point  (refer 
to  Figure  4).  Therefore,  during  the  process  of  solving  for  the 
the  boundary  layer  problem,  a  sliding  value  problem  on  the 
boundary  was  proposed.  That  is  it  is  possible  to  obtain  the 
sliding  value  (or  $n+2)  based  on  the  distribution  curve  of 
on  various  boundaries.  The  actual  methods  can  be  obtained  by 
referring  to  Reference  [2]. 
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Because  in  the  calculation  of  the  sliding  value  we  must 
consider  the  distribution  of  u  on  the  boundary,  therefore  the 
shape  of  the  distribution  curve  of  u  on  the  boundary  would 
directly  influence  the  accuracy  of  the  results  of  computation 
of  the  boundary  layer  flow.  Because  the  u  value  near  the 
boundary  is  very  small,  from  the  set  of  conservation  equations 
(2),  (3),  and  (4)  we  can  see  that  the  convection  along  the  x  - 
axis  can  be  neglected.  This  transforms  and  simplifies  the  flow 
near  the  boundary  into  the  model  of  the  Couette  flow.  Based  on 
the  characteristics  of  the  Couette  flow  on  the  boundary,  we 
can  derive  the  value  of  u  and  the  exponent  of  the  distribution 
curve  of  the  value  of  4>. 


Figure  4.  Definition  of  Sliding  Value. 

Key:  I  boundary,"  2.  True  variation  of  $  value,  3.  True 

variation  of  <&  value,  4.  E  boundary. 

Regarding  the  treatment  of  the  source  term  -£•  :  For 

dx 

the  flow  with  given  iL  distribution,  we  can  directly  use 

dx 

the  above  method  to  solve  the  problem.  But  for  those  problems 
with  the  distribution  of  yet  to  be  determined  in  the  flow 

dx 

process  (such  as  the  boundary  layer  flow  problem  in  the  restric¬ 
ted  duct),  it  is  relatively  more  complicated.  Rigorously 
speaking,  the  value  of  the  distribution  of  should  be 


resolved  based  on  a  Iterative  computation  method.  This  requires 
to  save  all  the  parameters  relevant  to  it.  in  the  flow  field 

dx 

and  to  carry  out  multiple  iterations.  This  would  significantly 
increase  the  computation  time  and  memory  units  required  which 
eliminates  the  advantages  short  time  and  few  memory  units 
necessary  for  the  original  progressing  method  used  to  obtain 
the  solution  to  the  parabolic  equation.  For  that  an  approxi¬ 
mation  method  is  used. 

The  expression  for  ^  can  be  derived  based  on  the 
continuity  equation  and  momentum  equation  of  the  one-dimensional 
flow  as  follows: 


therefore 


where 


dp  __  _ jF|  _  2u  dih  ^  mu  f  _1_  dA  _1_  dp  \ 

dx  A  A  dx  A  'A  dx  p  dx  1 


dp  \  dp  dJT 

fi  9  ~  T 

F"  25  dm  diu  dA  rhu_  df 

dp  _  A  A  dx  +  A*  dx  A?  dx 

dx  j  _  w5 

AP 

u  -  J  rpu’dy/^  rpudy,  T  -  j  rpuTdy  /j  rpudy 


(23) 


(24) 


For  gas  flow  of  low  Mach  number,  the  effect  of  pressure 
on  density  can  be  neglected.  Thus  the  second  term  of  the 
denominator  on  the  right  side  of  the  equal  sign  no  longer  exists. 

As  for  the  physical  meaning  of  ^  ~  ,  it  can  be  looked 

at  as  a  discussion  of  a  compressor  flow  problem.  When  the 
computation  progressed  to  x  ~  x.  ,  it  was  found  that  the 
results  were  that  the  fluid  did  not  fill  the  entire  cross-section 
of  the  compressor.  This  indicated  that  the  original  given  it 
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value  was  incorrect.  We  should  have  started  with  a  new  assump¬ 
tion  of  4^  to  obtain  the  cross-section  of  the  fluid  until  it 
finally  matched  the  cross-section  of  the  compressor.  At  the 
present  moment  this  method  is  not  used.  Instead  the  computation 
proceeds  to  x  —  xD  .  But  the  Ad  value  which  has  an  effect 

dx 

on  the  At  value  is  calculated  using  the  following  equation 

dx 

dA _ C Ad.  a  A,.  ,)  ( 25  ) 

dx  XD  —  Xm 

This  is  to  say  that  the  difference  in  the  At  value  at  the 

dx 

previous  stop  (x  —  x.)  is  partially  compensated  in  the  effect 
of  area  variation  at  the  next  stop  .  In  the  actual 

calculation,  a  partial  value  of  equation  (25)  should  be  used. 
Only  by  doing  so  that  the  instability  of  the  calculated  can 
be  prevented. 

In  using  the  above  method,  it  should  be  noticed  that  Af  D 
and  Ad  D  must  be  very  close  at  all  times.  They  should  be  the 
same  at  some  boundary  points.  Otherwise  the  approximation  of 
the  calculated  flow  field  has  no  meaning  whatsoever.  (See 
Figure  5). 
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Figure  5.  The  physical  meaning  of  dA/dx 

As  for  the  amount  of  partial  value  of  equation  (25)  that 
should  be  used,  it  depends  on  the  closeness  of  A^.  ^  and  A^  ^ 
which  is  checked  during  the  computation  process.  For  an  approxi¬ 
mate  straight  tube  flow,  we  can  use  0.1  -  0.2  as  the  correction 
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coefficient.  For  compressors  which  compress  quickly,  the 
correction  coefficient  is  very  small  or  is  given  by  sections. 


III.  The  Contents  of  the  Computation  Program  and 
Computed  Example 


Computation  program  has  been  compiled  based  on  the  above 
discussion.  Due  to  the  limitation  in  pages,  the  block  diagram 
and  the  actual  content  of  the  relevant  program  are  omitted.  The 
computer  used  to  process  the  program  was  a  Felix  C-256  computer 


We  have  used  this  program  to  calculate  a  ring  shaped  com¬ 
pressor  (Figure  6).  With  regard  to  the  axial  static  pressure 
variation,  the  calculated  and  experimental  values  were  compared 
and  shown  in  the  figure.  The  results  are  pretty  close. 
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Figure  6.  Computer  Example.  Key:  1.  p/p  inlet,  2.  L/L 
total  length,  3.  direction  of  flow,  U.  symmetry  axis,  5.  experi 
mental  value,  6.  calculated  value,  7.  (a)  axial  static  pressure 
variation  of  the  compressor,  9.  (b)  geometric  diagram  of  the 
ring  shaped  compressor. 


IV.  Closing  Remarks 


This  paper  based  on  the  momentum  equation  expressed  by 
stress  tensor  published  in  Reference  [6]  to  derive  the  momen¬ 
tum  equation  of  non-rotational  steady  axi-symmetric  flow  of  a 
compressible  fluid  in  the  turbulent  boundary  layer  in  the  axi- 
symmetrical  coordinate  system.  It  can  be  found  very  clearly 
from  the  derivation  and  discussion  that  equation  (2')  (similarly 
equations  (3')*  (V)  and  (18))  is  not  only  suitable  in  solving 
the  boundary  xeyer  problem  on  the  surface  of  a  revolved  body 
but  also  appl.  aMe  to  the  boundary  layer  problem  of  an  axi- 
symmetric  jet  flow. 

On  the  basis  of  Reference  [2],  the  computer  program  after 
being  partially  modified  by  us  can  be  applied  to  the  following 
two  situations  from  actual  computational  verification. 

1.  When  the  entire  flow  in  the  duct  is  treated  by  taking 
the  viscosity  into  consideration,  if  the  curvatore  of  the  wall 
surface  of  the  duct  is  not  too  large  and  the  rate  of  curvature 
change  in  the  direction  of  the  flow  is  also  not  too  large,  then 
the  entire  flow  field  can  be  calculated  based  on  this  method. 

2.  If  the  curvature  change  of  the  duct  is  large  then  it 
is  necessary  to  combine  this  method  and  the  non-viscous  flow 
in  the  main  flow  method  to  compute  the  solution. 

Appendix  I 

The  Deviation  of  the  Continuity  Equation  and  the  Momentum 

Equation 

When  the  compressible  viscous  fluid  flows  steadily  in 
space,  its  continuity  equation  and  momentum  equation  (without 
external  force)  are  as  follows: 


v  •  0>V)  -  0 


(1-1) 


▼  (y  |VI*)- YX  (v  Xv)-^v«  (1-2) 

Based  on  the  vector  and  tensor  operations  introduced  in 
Reference  [6],  for  the  non-rotational  axi-symmetric  flow  in  the 
axi-symmetrical  coordinate  system  (Figure  1)  considering  that 
ks  and  are  very  small,  then  *. -i,  *.-i,  and  . 

The  coordinates  corresponding  to  m.e.r  are  *,  y,  e  .  The  com¬ 
ponents  correspond  to  the  vector  quantity  V  are  and 

»,-o  .  In  addition,  all  the  variations  in  the  9  direction 

for  all  the  variables  can  be  neglected1- 

Therefore,  the  continuity  equation  can  be  written  as  the 
following: 


(1) 


The  equation  of  the  x-direction  momentum  component  is 


du 


.  a« 


Ar.,  a*,,  tin 3  cot  9  tin 3 

+  73?  +  —  T  +  ***-7“  ~  — 


The  relation  between  the  stress  tension  it  and  the  viscous 
stress  tensor  t  is  as  follows: 

*11  m  *ll  +  *Ht 


where  61j  is  the  Kronecker  5,  we  get 

du  * «  I  dp  I  *dr„  i  T.,  .  T.,  tind  U,  cm  3  :»  tin  3 

"a*  +  ,'5ym~Pd*+P  a*  +  P  dy  +  p  r  P  C  -  ' 


Based  on  the  qssumption  made  by  Prandtl  on  the  boundary 
layer,  the  above  equation  can  be  simplified  into 


du  du 


dy 


cot  d  dp 

~  *■»  -  £ 


(1-3) 


19 


After  considering  the  effect  of  curvature,  we  can  obtain 
the  momentum  component  equation  in  the  y  equation  as 


tt  I 

Ku  “  7  a7 


Therefore 


Since  the  boundary  layer  belongs  to  the  o(s)  class,  the 
static  pressure  difference  in  the  y  direction  in  the  actual 
boundary  layer  can  be  ignored. 

Based  on  the  conditions  of  the  non-rotational  axi- 
symmetrical  flow  and  from  an  order  of  magnitude  comparison, 
we  know  that  the  momentum  component  equation  in  the  0  direction 
no  longer  exists. 


For  simplification,  let  t  represent  x  and  also  consider 

...  xy 


-*-T  and  ir--3T 


(because  the  static  pressure  difference 


in  the  y  direction  can  be  ignored).  Equation  (1-3)  can  then 
be  transformed  into  the  following: 


dii  dr  cot  9 
dy  “  dy  +  * 


(1-4) 


As  for  the  mass  transfer  equation  (3)  and  energy  equation 
(4),  they  can  be  derived  using  the  same  procedures  as  described 
above.  It  will  hot  be  repeated  here. 
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Appendix  II.  The  Principles  Used  in  the 
Derivation  of  the  Difference  Equation 

When  establishing  the  difference  equation,  the  relation 
between  and  *-»  used  were  linear  and  stair-case  shaped 

distributions,  respectively.  As  for  the  $  value  used  in  A. 

out 

the  $  value  at  the  downstream  xn  was  used.  This  could  insure 

u  r  o  1 

the  stability  of  the  difference  computation1-  J .  But  the  coeffi¬ 
cient  a,  b,  and  c  were  obtained  from  the  $  value  at  xy  upstream. 
The  difference  format  of  the  equation  was  then  established  based 
on  the  integration  with  respect  to  the  entire  control  volume  for 
all  the  terms  in  equation  (18). 

The  source  term  d  has  different  content  corresponding  to 
different  conservation  equations.  For  mass  transfer  and  energy 
conservation  equations,  we  assumed  that  in  the  entire  control 
volume  the  value  of  d  was  uniform  and  equal  to  the  value  dp 
(i.e.  d-uj  or  d-x  varies  as  a  stair  case  shaped  variation). 

For  the  momentum  conservation  equation,  the  corresponding 
value  of  $  is  u.  The  calculated  result  of  u  would  simultaneously 
influence  the  solutions  to  the  mass  transfer  and  the  energy  con¬ 
servation  equations.  Therefore  the  importance  of  the  more 
accurate  assumption  of  the  source  term  distribution  in  the  momen¬ 
tum  conservation  equation  becomes  more  apparent.  Let  us  assume 
that  the  variation  of  d-w  between  network  nodal  points  was  linear 
and  the  variation  of  d-x  between  the  network  nodal  points  was  a 
staircase  shaped  variation. 

Based  on  the  above  principles,  we  can  derive  the  difference 
equation  of  the  mass  transfer  or  energy  conservation  equation 

l,»o*  +  *••(>  +  l»#u-  *  t>m  t.(*o+  -  «D)  -  fl(«0  ~  ♦»-)  rf.  +  ("fs ^,(*0  ~  *•)  (  19  ) 
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where 


1  (lap*  -  <Qp)  •  .  *  0*o+  +  3tt>o) 

**  m  *4  (*o  —  ”  *>0")  —  Wo*)  4  (“o*  -  «0*) 

3  4_ 

*'  ”  4^xfl  —  *.)  *"  4 

I  (a>p  -  a>D-)  «  _  (aip-  4-  3a>p) 

•»  “  7  (*o  -  ».X«o+  -  »o')  ”  ("?♦  ~  «o-)  ~  T  (»o+  -  "o-) 

*  (a>n»  -  a>o)  ..  .  3 

*•"“4  (*o  -  -«o-)  *•*  ”  4(x0  -  ».)  * 

(<»0  -  *0-)  J 

4  (xo -x.X»o4-«o-J  ' 

2  C„+ _ 

—  ®>0“XWO+  ”  Wb^ 

2C..- 

*'  "*  (tfo+-Wo-)(wO- 


Similarly,  the  difference  equation  of  the  corresponding 
momentum  conservation  equation  can  be  derived 


*, uo*  +  *i«o  +  tt*o-  +  i,  “  8tC*o+  -  “a)  ~  tt(" o  ~  *o~)  +  *.«>>♦  +  3«“o  +  i*“°"  +  5* 


where  e„  t„  e>'.  *.  *.  s»  were  defined  as  above,  and 
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Abstract 

A  momentum  equation  of  .steady  two-dimensional  axi-symmetrical  turbulent  bound¬ 
ary  layer  is  derived  employing  stress  tensor  analysis.  It  is  found  that  if  the  flow  in 
a  duct  is  to  be  treated  as  viscous  everywhere,  the  equation  is  valid  only  if  the  curvature 
of  the  duct  and  the  rate  of  change  of  the  curvature  of  the  duct  in  the  flow  direction 
are  small.  These  requirements  so  far  have  not  been  verified  by  previous  authors. 

In  this  paper,  a  computer  program  for  solving  the  governing  equations  has  been 
completed  and  a  numerical  example  is  selected  to  show  its  degree  of  accuracy. 
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LAMINAR  HEAT  TRANSFER  WITH  MASS 
INJECTION  AND  CHEMICAL  REACTION 


Zu  Tie-lin1 


In  order  to  simplify  the  boundary  layer 
problem  with  mass  injection  and  chemical  reaction, 
a  general  stoichiometric  formula  has  been  derived. 
Using  the  chemical  equilibrium  as  an  example  we 
performed  an  analysis.  The  calculated  results  were 
found  to  be  consistent  with  the  experimental  ones. 
When  the  amount  of  injection  is  zero,  the  relevant 
data  were  in  good  agreement  with  the  results  listed 
in  Reference  [4],  Finally,  it  was  pointed  out  that 
with  the  increasing  available  energy  of  the  exo¬ 
thermic  reaction  the  effect  of  Lewis  number  Le  on 
the  heat  transfer  decreases. 


In  the  boundary  layer  problem  with  mass  injection  and 
chemical  reaction,  one  usually  assumes  that  the  Lewi  ^  number 
Le  =  1  in  the  combustion  loss  calculation  in  order  ^  o  obtain 
relatively  simple  results.  Along  with  increasing  amount  of 
injection  and  variation  of  other  relevant  conditions,  this 
assumption  would  cause  significant  errors  in  the  results  of 
calculated  heat  transfer  and  effective  combustion  heat.  Lees^-1^ 
has  provided  an  analytical  solution  for  Le  ¥  1  which  is  expressed 
by  Blasius  function  to  the  wall  surface  chemical  equilibrium 
problem  of  a  frozen  boundary  layer  with  a  constant  transfer 
characteristic.  This  paper  studied  the  affect  of  chemical 
reaction  on  the  heat  transfer  through  a  discussion  of  a  chemical 
equilibrium  boundary  layer  with  Le  ¥  1.  We  have  carried  out  a 
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discussion  on  the  equations  of  boundary  layer  chemical  reaction 
and  wall  surface  conditions  and  obtained  an  universal  generali¬ 
zed  stoichiometric  formula.  Using  the  assumption  imposed  by 
Lees,  we  provided  the  numerical  solution  to  a  series  of  examples 
for  the  equilibrium  boundary  layer.  After  further  making  a 
linearity  assumption,  we  obtained  an  approximate  equilibrium 
solution  similar  to  Lees'  frozen  solution. 

In  the  constant  transfer  assumption,  the  most  worthwhile 
discussing  subject  is  the  applicability  of  using  Fisk's  Law  to 
express  the  mass  flow  and  taking  Le  as  a  constant.  The  author 
has  written  an  article  on  this  subject  with  Yao  Kang-Chun  and 
Hu  Cheng-Hwa  in  1964  (unpublished).  That  paper  compared  the 
heat  transfer  calculated  in  the  decomposition  of  air  using  the 
multiple  element  method  and  the  two-element  method  with  a  con¬ 
stant  Le  number.  Its  conclusion  was  that  when  the  constant  Le 
number  was  properly  chose.)  we  could  obtain  the  same  heat  trans¬ 
fer  as  the  one  calculated  by  the  multiple  element  method. 

(Figure  1).  As  for  the  combustion  reactions  in  decomposed  air, 
as  long  as  a  proper  Le  number  is  chosen,  we  will  get  the  same 
result . 


I.  Basic  Equations  and  Generalized 
Stoichiometric  Formula 

The  continuity,  diffusion,  momentum  and  energy  (expressed  by 
the  frozen  enthalpy  HT Equations  of  a  steady  laminar  boundary 
layer  with  chemical  reaction  are: 

(?«»■•),  +  (prr*),  —  0 
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Figure  1.  Heat  Transfer  Calculated  Using  the  Multiple 
Element  Method  and  Two-Element  Method.  Key:  2.  Multiple 

element  method,  3.  Two-element  method.  _ — -jgjc 
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Where  r  is  the  radius  of  the  revoluted  body;  when  e  *  0 
it  is  a  two  dimensional  flow;  When  e  ■  1,  it  is  a  three 
dimensional  axi-symmetrical  flow;  p  is  the  density  of  the  mass 
u  and  v  are  respectively  the  velocity  components  in  the  x  and 
y  direction;  x  and  y  are  respectively  coordinates  along  the 
direction  of  the  object  and  that  perpendicular  to  the  surface 
of  the  object;  the  subscript  i  represents  any  element  in  the 
complete  element  1;  K  is  the  concentration  (mass  ratio);  J  is 
the  mass  flow;  w  is  the  chemical  rate  of  formation;  P  is  the 


m 

pressure;  u  is  the  viscosity;  the  frozen  enthalpy 

k  »  K,ht,  h,  ■»  I  ct,dT\ctt  is  the  isobaric  specific  heat  of  element 

t 

i;  the  isobaric  specific  heat  of  the  gas  mixture  > 

K  is  the  thermal  conductivity;  the  Praneltl  number  Pr  *  £*£  ; 

0  * 

h,  is  the  enthalpy  of  formation  at  0°K  for  a  unit  mass  of 

element  i.  _ _ 

Let  us  assume  that  every  element  is  a  perfect  gas,  then 

FMi-etRT  (5) 

PiZ-pRT  (6) 

When  M  «  K,  »  and  is  the  molecular  weight  of  the  ele¬ 
ment  i;  R  is  the  conventional  gas  constant.  Let  us  assume  that 
the  mass  flow  can  be  expressed  using  Fisk's  Law,  which  treats 
the  diffusion  coefficients  D^j  of  the  two  elements  the  same 
(both  to  be  D)  and  only  considers  concentration  diffusion,  then 
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When  is  frozen  at  zero  and  it  is  a  function  of  the  local 
p,  T,  and  K.^  when  there  is  chemical  reaction.  The  functional 
relation  is  given  by  the  chemical  kinetic  conditions.  Let  us 
assume  that  there  are  a  total  number  of  s  independent  chemical 
reactions  in  the  problem  we  are  considering: 


jXx-o  (8) 

<•1 

where  is  the  stiochiometr ic  coefficient  of  element  i  in  the 
pth  chemical  reaction;  is  the  molecular  formula  of  element 
i.  Every  reaction  in  equation  (8)  has  only  one  independent 
reaction  rate,  s  is  usually  smaller  than  1.  In  order  to  solve 
for  the  mass  ratio  of  the  entire  1  elements,  the  1  -  s  diffusion 
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equations  necessary  to  be  preserved  are  usually  transformed 
into  the  frozen  form  using  the  conservation  of  element  method. 
This  method  is  mandatory  for  equilibrium  problems.  For  non¬ 
equilibrium  problems,  the  calculation  procedure  can  be  simpli¬ 
fied  and  it  is  irrelevant  whether  the  mass  flow  is  expressed 
using  equation  (7).  In  order  to  overcome  the  disadvantages  of 
using  the  conservation  of  element  method  such  as  the  inconven¬ 
ience  to  put  it  into  a  routine  form  and  inability  to  choose 
the  reference  element,  this  paper  obtained  a  "generalized 
stoichiometric  formula." 


Let  us  make  the  reaction  rates  of  1  -  s  elements  out  of 
the  1  elements  in  the  1  -  s  reactions  »},»],••♦,»;  to  be  indep¬ 
endent,  then  the  rate  of  formation  of  each  element  can  directly 
be  expressed  oy  the  linear  combination  of  the  s  independent  rates 
of  formation: 

(9) 

where 
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If  we  consider  the  problem  as  a  linear  space  problem,  then  this 
set  of  independent  rates  of  formation  can  be  mathematically 
treated  as  a  set  of  independent  vectors  of  a  "basis".  The 
coefficient  a^  then  becomes  the  coordinate  under  this  "basis". 
The  appearing  in  equation  (9)  is  not  convenient  to  use,  it 
is  necessary  to  select  a  new  "basis"  («»,  »»,*••»  «*.)»  which  is 
expressed  by  the  total  reaction  rate.  The  reaction  rate  not 
included  in  the  new  "basis"  or  the  non-independent  reaction 
rate  can  be  expressed  linearly  using  the  new  "basis": 

* 

«*«  “  2  *1“/  '  +  2,  •••,  /) 

<•#•1 


(11) 


where 


(12) 


Equation  (11)  is  the  generalized  stoichiometric  formula.  The 
use  of  this  formula  is  especially  convenient  for  problems 
with  complicated  chemical  reactions.  Let  us  define  the  symbol: 

s-+';i-^(c'i)  <13> 

and  then  equation  (2)  can  be  written  as: 


Substituting  equation  (14)  in  equation  (11),  we  get 

•  0  (15) 

2*1*1  <*-'+  !,•••, O  (16) 

<-! 

Equation  (15)  is  the  preserved  (1  -  s)  diffusion  equations  in 
the  frozen  form.  The  remaining  s  equations  are  provided  by 
the  s  chemical  reaction  conditions.  When  the  chemical  reactions 
are  in  equilibrium,  the  s  conditions  are 

(17) 

4-1 

f*  V*i 

where  JX',  is  the  balancing  constant  of  the  p  reaction.  If 
we  slightly  change  the  operating  symbol  by  substituting 
with  the  rate  of  formation  per  unit  area  m^,  then  equation  (15) 
can  be  expanded  to  the  boundary  or  any  plane. 


II.  Approximate  Solution  and  Iteration 
Equations 

Let  us  consider  pr,  Le  (  -  pDe,/0,  /(  -  pmAwO  to  be  reference 
constants.  For  practical  reasonsj  let  us  limit  our  discussion 
to  the  point  of  the  revoluted  body.  The  following  type  of 
correction  was  made  using  the  Mangler-Dorodonitsyn  transforma¬ 
tion: 


(18) 


2 Vs  * 

e  ~  J#pww,4f* 

where  the  subscripts  e  and  w  represent  the  wall  surface  and 
outer  fringe  of  the  boundary  layer  values,  respectively.  Let 
us  define  the  following  dimensionless  quantities: 


'■  ■ 3  *  ~  ~  l  "A:  ~‘-dT- 


For  cold  walls,  we  can  neglect  the  pressure 
In  the  momentum  equation.  Let  us  assume  that 
element  Is  the  same  and  use  equation  (11)  in  the 
tion,  then  the  basic  equations  can  be  written  as 


gradient  term 
the  0^^  for  every 
energy  equa- 


■  0 

•  (19) 

+  Sc -■  o  -»/  +  l, /  +  /> 

(20) 

(tr  +  Ufl  +  W  +  L«  2 

-  Pr/CLe  -  1)  5]  SUK,. 

(21) 
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* 

boundary 

where 


conditions  for  these  equations  are: 


/,<o)  -  o,  m  -  r-rU.i*) 

/,(«)- 2  .  .  . 

**(0)  - 

/r(®)  "  |r»»  /xC00)  "  1 


(22) 


As  for  the  s  equations  lacking  in  finding  solutions  to  K^, 
equation  (17)  will  be  used  as  the  supplement  for  equilibrium 
problems. 


We  carried  out  integration  for  equations  (20)  and  (21)  and 
wrote  the  result  of  integration  of  equation  (21)  in  an  iterative 
form.  We  get 


^  (23) 

“  (£*.  ^*»)^(®i/»»®c)  ( 24  ) 

<*7  +  Lct»y>  -  (,T  +  LeJ.0.  -  <ir.  -  It.  +  [Le  -  (Le  -  1)  G'-'>(oo)]  ( 25  ) 

•  a,w  -  /.,  Pr)  +  G(m~‘Kv) 

igi  +  Le^cOw  “  {gn  —  gjm  +  [Le  —  (L«  —  1)G<*“',(00)]C^ «»»  —  ^-w)}  (  26  ) 

•  £,(0;  /„,  Pr) 

where  B  is  the  Blasius  function: 

9  OjJ  f~, «)  -  J*  IWn /(,  fWn 

z  =  Sc  or  Pr.  The  detailed  numerical  table  of  B  and 
given  in  Reference  [2]. 
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Under  the  situation  that  the  boundary  conditions  are  known, 
equation  (26)  can  provide  the  heat  transfer  which  we  are 
interested  in  obtaining.  The  G(»)inthe  equation  usually 
has  to  be  obtained  through  the  iteration  equation  set  (23)  - 
(29)  and  equation  (17).  The  first  order  approximation  of  G(») 
can  be  calculated  using  the  when  Le  =  1.  Calculation  makes 
clear  that  convegence  in  this  form  of  equation  is  very  quick. 

When  the  temperature  in  the  boundary  layer  exceeds  2000°K, 
because  the  vibrational  degrees  of  freedom  for  every  element 
are  nearly  all  excited,  Cp  can  be  treated  approximately  as  a 
constant.  At  this  time  in  equations  (25)  and  (26)  — 0  — jyr,  ■ 

It  is  even  simpler  to  find  the  solution. 

III.  Wall  Surface  Mass  Ratio 
and  Heat  Transfer 

By  integrating  equation  (20)  from  n  =  0  to  jj  —  o  +  *(• 
is  a  small  quantity)  and  also  noticing  equation  (24),  we  can 
obtain  the  conservation  of  mass  conditions  of  wall  surface 

-  (£A— a  +  £*.)/(  s  +  0  ( 30 ) 

where  j?4  is  defined  by  equation  (16);  Klw  -  represents  the 
mass  ratio  before  injection  (i  =  K  or  j),  for  non-injection 
elements  Klw_  *  0;  and  Kiw  represents  the  mass  ratio  after  the 
injection  and  chemical  reaction. 

B~~  (31) 

Using  wall  temperature  as  a  parameter,  equation  (30)  and 
equation  (17)  (or  general  catalytic  conditions)  form  a  set  of 
closed  equations  within  which  K^w  can  be  solved. 

With  different  mass  injection  conditions,  we  can 
divide  the  wall  surface  onto  three  types  of  problems  and  they 


are  the  "homogeneous  reaction",  the  "heterogeneous  reaction", 

and  the  combination  of  the  two  above.  For  "heterogeneous 

reaction"  (such  as  the  combustion  of  carbon),  the  quantity  B 

of  interest  in  the  combustion  loss  calculation  can  be  directly 

solved  using  the  above  method.  This  is  because  the  K,  of  the 

iw 

solid  phase  element  is  always  1  in  equation  (17)  and  it  is  0 
in  equation  (30-  For  "homogeneous  reaction",  the  mass  ratio 
before  injection  should  be  given  by  the  condition  in  the  wall 
(material  composition,  control  method,  etc.). 

Using  the  approximation  method  to  integrate  the  energy 
equation  (21)  near  the  wall  surface,  we  can  get  the  heat  trans¬ 
fer  into  the  wall  as 


Pr)/*Pr-‘ 

Oy  _ 

x  ( w* )«  (h"  ~  htw 

+  L*rAAth  -  B  LeriLE~) 


where 


AA,ii  -*  2  {Ki, —  [(B  +  l)K/„  (33) 

W  -  1  +  <Lc  -  1)10  -  G(«5))  +  «,(C( oo)  -  G,(°o»]  (  34  ) 

y  *»  1  —  G(oo)  +  H*[G( oo)  —  G/(oo)]  (  35  ) 

The  latent  heat  of  phase  change  is  L£  *  HTW  -  HTW_.  The  sub¬ 
script  f  represents  the  frozen  value.  When  frozen,  c(oo)  —  G,(oo), 
Ler/,  n  *»  1  —  G,(oo)  .  was  given  in  Reference  [1]. 

The  subscript  s  represents  the  stationary  point  value. 


If  we  neglect  the  differences  of  parameters  such  as  Pr,  Le , 
1,  etc.  between  the  equilibrium  and  frozen  conditions,  from 
equation  (32)  we  can  see  that  the  exponent  y  is  the  only  para- 
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Igure  2.  Comparison  of  Theoretical  and  Experimental  Results 
Key:  2.  Experimental  points  (arranged  based  on  Le  =  lC5]), 

3.  Experimental  points  (arranged  based  on  Le  -  1.4),  4.  Theo¬ 
retical  curve  (Le  =  1),  5.  Theoretical  curve  derived  in  this 
work  (Le  =  1.4),  6.  Heat  transfer  of  carbon  comhustion  in 

decomposed  air,  7.  Heat  transfer  of  decomposed  air  under  same 
boundary  conditions. 


meter  to  indicate  the  effect  of  chemical  kinetic  condition  in 
the  boundary  layer  on  the  heat  transfer.  Similar  to  the  kinetic 
energy  recovery  factor  which  depends  on  the  Pr  number,  Ler  here 
corresponds  to  the  extent  of  the  transformation  of  the  excess 
chemical  enthalpy  of  the  outer  fringe  compared  to  the  wall  sur¬ 
face  into  thermal  enthalpy.  .  Because  the  generalized  stiochio- 
metric  formula  used  can  arbitrarily  choose  the  reference  ele¬ 
ment,  it  is  very  convenient  to  select  the  element  with  concen¬ 
tration  at  the  wall  surface  to  be  zero  (not  zero  in  the  outer 
fringe)  as  the  reference  element  in  the  calculation  of  equations 
(32)  -  (35). 
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e1fnd6fhP  distributions  of  the  dimensionless  temperature 

and  the  nitrogen  and  oxygen  atomic  mass  ratios  along  n  in 
the  equilibrium  boundary  layer  of  decomposed  air.  S 

^this’work1  SUrfaCe  reacti°n  equilibrium,  3.  this  work. 
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figure  4.  The  distributions  of  the  dimensionless  temperature  9 
and  the  mass  ratio  k.o-c o„  co,o>  along  n  in  the  combustion  of 
carbon  monoxide,  r. .  6oook  t.  -  2oook’p  -  3.5am,  -/. -0.25  Le  -  j.2  pr  -  0.71 
Key:  2.  wall  surface  reaction,  3.  equilibrium:  the  reaction  in 


the  boundary  layer  is  CO  +  0  =  C02  (assuming  nitrogen^ does  not' 


recombine),  4.  equilibrium  distribution,  5.  frozen  distribution. 


This  paper  conducted  heat  transfer  calculations  for  the 
following  examples.  1)  The  homogeneous  injection  of  air,  the 
reaction  equation  is  0,«-fc20s  N2=*—  2N  ;  2)  the  homogeneous 

injection  of  oxygen,  which  is  the  injection  of  oxygen  in  an 
scattered  flow  of  oxygen,  the  reaction  equation  is  0,«^-*20  ; 

3)  the  homogeneous  injection  of  nitrogen,  the  reaction  equation 
is  N,«*-*2N  ;  4)  the  injection  of  carbon  monoxide,  the 

reaction  equation  is  co  +  O^^COj  »  5)  the  injection  of  carbon 
monoxide,  the  reaction  equations  are  CO  +  O^^CO,,  20^—0,  ; 

6)  the  injection  of  carbon  monoxide,  the  reaction  equations  are 
-w^CO,,  20*-*Oj,  2N<*-*N,  •  The  results  of  example  1  at  zero 

injection  agreed  with  the  results  of  Reference  [4]  very  well. 

In  Reference  [5],  they  have  found  that  the  theoretical  results 
obtained  by  taking  Le  =  1  was  far  different  from  the  results 
they  obtained  from  their  carbon  combustion  experiment.  If  the 
original  data  in  Reference  [5]  were  rearranged  using  the  method 
described  in  this  paper  and  let  Le  =  1.4,  then  it  was  found 
that  all  the  experimental  points  were  located  either  above  or 
below  the  theoretical  values  obtained  using  this  method  (see 
Figure  2).  Figures  3-6  provided  the  distributions  of 
v  ~  __  v'  <£,  t,  fJ  du  \‘  along  n .  Figures  7-8  showed  the  varia¬ 


tions  of  fw,  Tg,  and  Tw  with  y.  The  two  extrema  in  the  distri¬ 
bution  of  •y'uixi  A  / du,\  along  n  correspond  to  the  maximum 

/  m  |  P  •  '  dx  •  * 


exothermic  and  endothermic  surfaces. 


T.  -  6000K  r.  —  1J00K  P-3.5»ln» 
o  #|1  A  *2  v  «|3  Off* 

H  r  <•  '* 


Figure  7.  The  variation  of  the  chemical  kinetic  parameter 
y  with  -  V 

Key:  1.  frozen  value,  2.  example  4,  4.  example  1,  5-  example  2 

6.  example  3,  7.  example  4,  8.  example  1,  9-  example  2, 

10.  example  3. 


Figure  8.  The  variation  of  the  chemical  kinetic  parameter  y 


with  the  outer 
Key:  4.  example 


fringe  temperature  T  . 

C-/.-0.3) 

- r.  -  1500K1 .  ,  . 
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O7.-  1500K 


IV.  The  Approximate  Solution  and  Discussion 


When  Le  is  not  too  far  away  from  1,  we  can  expand  and 
gT  around  (Le  -  1)  and  obtain: 

K,  -KP+  (Le  -  DKJ*>  +  0[(Lc  -  J)»]  (36) 

ft  -  /?>  +  (Le  -  l)f?»  +  0((U  -  O’] 

(37) 

Let  us  neglect  all  the  small  terms  after  the  second  order  term 
of  (Le  -  1),  then  G(«)  can  be  solved  from  the  energy  equation 
with  Le  *  1.  From  equation  (21)  it  was  found  that  the  assump¬ 
tions  that  all  C  .  are  identical  and  that  -71  —  1  (or  mf4 — 1  ) 

pl  M  tf 

and  Le  -  1  are  small  quantities  of  the  same  order  of  magnitude 
are  equivalent.  The  concentration  of  the  selected  reference 
element  at  the  wall  surface  is  zero.  When  »&>n>0  and  m  <  >1  <  *>» 
*  is  the  value  of  n  at  6  =  0.99),  £*,-*0  and  is 
monoclinic  in  the  region  .  At  this  time,  H,  -  0,  Le' 

has  the  same  expression  as  in  the  frozen  state: 

Lc'  -  l  +  (Le  -  00  ~  G) 


or 


7  1  —  G 


Let  us  rewrite  int  0  -  KitSaB •  Where  S<  -  K,/Kh,  d  «.  (0  -  0,)/(0,  -  $t ), 

and  0^  and  02  are  the  dimensionless  temperatures  corresponding 
to  and  n2,  respectively.  Let  us  assume  that  St j  — 1  .  This 

assumption  is  suitable  for  the  air  decomposition  problem  when 
the  pressure  is  not  too  high.  For  other  combustion  reactions, 
it  corresponds  to  an  assumption  that  reaction  only  takes  place  at 
and  n2  (corresponding  to  the  maximum  exothermic  and  endo¬ 
thermic  surfaces  in  Figure  6).  Assuming  Pr  =  1  and  using  the 
subscripts  1  and  2  to  respectively  represent  the  values  at  n1 
and  n2>  we  obtain  the  following  by  integrating  Gj  between  the 


regions  to  n2: 


G'  "  !1d^J  +  7  w  -  /d]/C^  -  A.)  ( 38 ) 

Based  on  the  reaction  initiation  and  termination  conditions 
of  the  reference  element  at  and  n2,  the  n  value  is  determined 
using  the  method  of  determination  of  the  boundary  layer  thick¬ 
ness.  Furthermore,  through  Le  =  l  and  using  the  energy  equation 
with  constant  cp  as  well  as  equation  (17),  and  Bj2  can  be 

solved.  From  here  we  can  locate  the  corresponding  f  and  f 

nn 

from  the  Blasius  table.  The  calculation  indicates  that  the 
error  is  less  than  5%  comparing  the  y  value  obtained  using  the 
approximation  method  used  in  this  section  with  the  calculated 
results  obtained  in  the  last  section. 

When  fa-+Q,  fa— l  ,  G/—C,  (Gj  is  consistent  with  the 

results  in  Reference  [1]).  When  and  -*  1  also, 

G,—  l,  and  ti  —  o  .  This  corresponds  to  the  overlapping  of 
two  reaction  surfaces  at  the  outer  fringe.  Similar  to  the 
turbulent  layer,  a  concentration  type  of  discontinuity  appears 
at  the  outer  fringe.  The  transfer  process  has  no  effect  on  the 
heat  transfer.  In  reality,  situation  can  only  exist  close  to 
that  condition.  For  example,  when  Tg  =  6000  K  (cold  wall  condi¬ 
tion)  the  of  nitrogen  is  very  close  to  the  outer  fringe.  It 
is  meaningful  to  have  identical  results  when  G^  *  1  and  Le  =  1. 
If  we  define  the  extent  of  exothermic  reaction  occurring  in  a 
relatively  high  temperature  region  as  the  reactivity,  then  the 
assumption  that  Le  =  1  is  only  suitable  for  problems  with  strong 
reactivity.  When  fa  —  0  and  fa  -»o  also,  we  can  get  G*  -• 

This  corresponds  to  the  situation  that  the  concentration  inter¬ 
ruption  surface  is  right  on  the  wall  surface.  When  the  injec¬ 
tion  quantity  is  constant,  the  y  value  is  maximum  at  this  time 
which  also  means  that  at  this  time  the  effect  of  the  transfer 
process  has  the  most  effect  on  the  heat  transfer.  When  the 

40 


pressure  is  lower  and  the  temperature  of  the  outer  fringe  is 
higher,  the  oxygen  decomposition  reaction  is  approximately  in 
this  situation.  Corresponding  to  the  -fw  value  of  0.25,  0.50, 
0.75  and  1.00,  the  value  of  Gi -/«/£,( 0;/„)  are  0,  -0.51,  -1.52, 

-4.0,  and  -14.1,  respectively;  and  the  corresponding  Yj  are 
1.53>  1.84,  2.34,  3-33,  and  6.55  times  the  frozen  value  Yjf, 
respectively.  This  explains  that  when  there  is  chemical  reaction 
going  on,  y  may  also  be  greater  than  Yf. 

When  0<— /„<0.5  using  0^  and  0^2  as  variables  we  can 
linearly  expand  G^  around  the  two  points  0^  =  0  and  0j2  =  1 
to  get 


c,  -  c,  + 


b-£]'« 


+  (1  ” ■  G/Xfti  O 


In  the  above  equation,  because  Gf  is  less  than  1  and  also  less 
than  | fw |  and  fw  is  always  negative  in  injection  problems,  the 
coefficients  of  8^  and  0^2  are  always  positive.  Therefore,  Gj 
increases  with  increasing  0j 1  and  8j 2  and  Yj  decreases  with 
increasing  0^  and  Bj2.  This  explains  that:  when  the  reacti¬ 
vity  is  stronger  or  the  reaction  surface  is  closer  to  the 
higher  temperature  outer  fringe,  the  effect  of  the  transfer  pro¬ 
cess  on  the  heat  transfer  is  weaker.  All  the  factors  which  can 
increase  the  reactivity  such  as  reactions  with  high  equilibrium 
constants,  increasing  the  injection  quantity  of  the  reactant, 
decreasing  temperature  for  a  fixed  reaction,  increasing  the 
pressure,  etc,  can  make  Yj  smaller  (refer  to  Figures  7-8). 

For  problems  containing  two  or  more  reactions,  we  should 
still  calculate  the  contribution  to  the  total  e.’thalpy  change 
from  the  enthalpy  difference  —  Ar,„)&k  for  each  reaction 

based  on  equation  (28).  The  higher  the  enthalpy  difference  the 
more  its  recovery  capability  contributes  to  the  total  recovery 
capability  LeY.  Since  the  mass  ratio  and  enthalpy  of  formation 


of  nitrogen  far  exceeds  those  of  oxygen  in  air,  the  y  of 
decomposed  air  with  chemical  reaction  is  usually  smaller  than 
Y-  when  the  outer  fringe  temperature  is  sufficiently  high. 

With  increasing  outer  fringe  temperature,  other  reactions 
taking  place  in  the  decomposed  air  would  further  promote  the 
deviation  of  y  in  the  less  than  y^.  direction.  The  combustion 
reactions  in  decomposed  air  are  merely  reactions  between  the 
combustible  material  with  the  relevant  elements  in  air.  The 
trend  of  variation  of  y  is  similar  to  that  of  the  above  discussed 
pure  decomposed  air  problem. 


V.  Conclusions 


Using  the  generalized  stoichiometric  formulafto  eliminate 
the  non-independent  reaction  rates  has  the  convenience  of 
arbitrary  selection  of  a  reference  element.  This  method  is 
more  routine  than  any  other  methods.  Therefore,  it  is  more 
suitable  for  problems  with  complicated  chemical  reactions. 


The  heat  transfer  can  be  calculated  using  the  set  of  given 
iteration  equations.  Calculation  and  analysis  both  indicate 
that:  the  increase  in  reactivity  is  the  movement  of  the  reaction 

to  the  outer  fringe  and  it  has  the  effect  of  decreasing  the 
influence  of  Le  number  on  the  heat  transfer.  For  decomposition 
of  air  and  combustion  reactions  in  decomposed  air,  the  effect  of 
Le  number  on  the  heat  transfer  is  less  than  that  at  the  frozen 
state. 

In  the  study  of  non-equilibrium  problems  located  between 
the  equilibrium  and  frozen  states,  the  difficulty  is  not  in 
solving  the  boundary  layer  equation  itself.  Rather,  it  is 
difficult  to  provide  the  accurate  reaction  rates.  When  the 
amount  of  injection  is  not  too  large  and  near  the  stationary 
point,  the  difference  between  the  equilibrium  and  the  frozen 
states  is  not  too  significant.  Therefore,  under  such  condition 


the  non-equilibrium  condition  does  not  have  to  be  considered. 
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Abstract 

To  simplify  our  problem,  a  stoichiometric  formula  is  derived  for  boundary  layers 
with  mass  injection  and  chemical  reaction.  Aa  an  example,  the  solution  of  chemical 
equilibrium  was  analized.  Results  of  calculation  were  found  to  be  in  close  agreement 
with  that  of  experiment  [5],  and  with  that  of  [4]  when  there  is  no  injection.  Finally, 
it  was  indicated  that  the  effect  of  Lewis  number  on  beat  transfer  decreases  as  the 
capability  of  heat  generating  reactions  increases. 


THE  MEASUREMENTS  OP  THE  STATIC  AND  DYNAMIC  STABILITY  DERIVATIVES 
OF  CONICAL  MODELS  IN  THE  SHOCK  TUNNEL 

Ma  Jia-huan,  Tang  Zhong-beng,  Zhang  Xiao-ping  and  Guo  Yan-ping 
(Mechanics  Research  Institute,  Academy  Sinica) 

Experimental  research  of  shock  tunnel  for  the  development  of 
hypersonic  gasdynamics  has  been  conducted  widely.  Due  to  its 
extremely  short  duration,  there  exist  certain  difficulties  in  the 
measuring  technique.  Hence,  it  is  limited  to  the  static  aspect 
when  the  stability  of  a  flying  vehicle  is  studied.  In  fact,  the 
oscillatory  motion  of  the  flying  vehicle  after  reentry  directly 
affects  the  aerodynamic  load  and  the  aerodynamic  heating.  Hence, 
the  study  of  dynamic  stability  is  a  very  essential  issue.  For  the 
experimental  investigation  of  the  dynamic  stability,  the  model  free 
flight  method  which  avoids  completely  the  disturbance  of  the  sting 
displays  its  unique  superiority.  Trial  efforts  were  made  in  this 
field  but  have  not  been  successful  due  to  the  short  duration  of 
the  shock  tunnel  [1],  To  extend  the  range  of  application  of  the 
shock  tunnel  on  the  one  hand,  and  to  initiate  the  experimental 
study  of  the  hypersonic  dynamic  stability  on  the  other,  we  spent 
some  effort  on  modeling  and  angular  measurement  and  obtained  1.5- 
2.0  cycles  of  pitching  angle  motion  8c  =  10°  and  11°  cones  in  a 
hypersonic  flow  of  *  9.0  with  the  model  free  flight  method. 
Through  data  processing  technique,  not  only  their  static  stability 
derivatives  are  obtained,  but  also  the  preliminary  result  of  the 
dynamic  stability. 

1 .  Experimental  facility,  measuring  technique  and  model 

The  experiment  is  conducted  in  a  JF-8  reflection  type  shock 
tunnel.  The  air  flow  in  such  a  tube  is  motivated  by  the  mixing 
and  combusting  of  hydrogen  and  oxygen.  The  Internal  diameter  of 
the  3hock  tunnel  is  150  mm  and  the  diameter  of  the  test  section  is 
1.2  m.  The  typical  operating  condition  for  the  experiment  is: 
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=  9.0,  Re^  =  1.6  x  10^,  and  a  quasi-steady  operating  duration 
of  about  10  millisecond. 


The  first  step  of  the  free  flight  method  for  force  measure¬ 
ment  is  to  hang  the  model  in  the  test  section  of  the  tunnel  with 
extremely  thin  nylon  wires  at  prescribed  initial  conditions.  As 
soon  as  the  initiating  shock  scans  through  the  test  section,  the 
nylon  wires  are  burnt  out  and  the  model  is  left  exposed  to  the 
aerodynamic  force  and  the  gravity  only,  without  any  other  support. 
The  free  flight  motion  is  then  satisfied.  The  motion  of  the  model 
at  this  instant  can  be  recorded  by  means  of  synchronized  high  speed 
photography.  By  analyzing  the  data  and  the  motion  of  the  model, 
the  aerodynamic  characteristics  can  be  determined.  The  "Strobokin" 
high  speed  flasher  with  flashing  frequency  of  f  =  5  kc/s  is  used. 
Each  single  flashing  pulse  is  about  1  sec.  The  motion  history  of 
the  model  is  recorded  by  the  revoling  drum  camera.  The  duration 
time  of  the  flash  is  controlled  by  a  timer  and  is  corresponding 
to  the  quasi-steady  operation  time  after  the  tunnel  is  started. 

In  order  to  identify  the  flow  state  corresponding  to  each  film, 
the  measured  signal  of  the  pitot  pressure  in  the  test  section  of 
the  wind  tunnel  and  flashing  signal  must  be  recorded  simultaneously 
on  the  oscilloscope  (see  Figure  1).  The  complete  set  up  of  the 
testing  system  is  shown  in  Figure  2. 


A  model  of  very  small  moment  of  inertia  of  rotation  is  one  of 
the  basic  requirements  to  obtain  relative  more  periods  of  angular 
motion  in  the  shock  tunnel  and  hence  to  obtain  the  dynamic  deriva¬ 
tive.  During  the  quasi-steady  operating  time  t  of  the  wind  tunnel, 
the  period  of  oscillation  of  the  model  is 


2*V 


(1) 


and  within  the  same  duration,  the  flight  distance  of  the  model  due 
to  drag  is 
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By  direct  substitution,  it  can  be  observed  that 


N-J.  l~C-.  — 

*  V  Co*  T  1 


<Ut£ttBJR«n 


Figure  1.  Pitot  pressure 
and  flashing  signals 
upper  line:  flashing  impulse 
signal 

lower  line:  pitot  pressure 
signal 

scanning  speed:  2  ms/cm 


Figure  2.  Schematic  diagram 
of  testing  system 

I - waiting-type  drum  camera;  2-photoelectric 
diode;  3-  preamplifier;  4-preamplifier; 

5-  two-channels  oscilloscope;  6-  pitot  head; 

7-  transducer;  8-  trigger  delayer;  9-  pre¬ 
setting  circuit;  10-  digital  freq.  meter; 

II- flash  light;  12-  source  box;  13-  controller; 
14-  timer;  15-  electric  source; 


In  the  above  equations. 

t  quasi-steady  operation  time  duration 

q  dynamic  pressure  of  the  flow 

d  characteristic  dimension  of  the  model,  that  is,  the 
base  diameter  2 

S  the  base  area  of  the  model,  S  =  77 r- 

w 

m  the  mass  of  the  model  ,  m  —  - 

I  the  rotational  moment  of  inertia  of  the  model  about 
the  lateral  axis  through  the  center  of  gravity 

C  the  static  stability  derivative  of  the  model 
ma  J 

Cn  --.the  effective  drag  coefficient  of  the  model  under  the 
test  condition 
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The  model  should  be  designed  such  that  r  does  not  overshoot 
the  range  of  the  camera  lens  while  maintaining  the  maximum  poss¬ 
ible  periods  of  angular  motion.  From  Equation  (3),  it  is  ob¬ 
served  that  for  a  certain  aerodynamic  profile  under  a  specified 
flow  condition,  the  model  mu&t  acquire  a  sufficient  amount  of 
mass  with  a  minimal  rotational  moment  of  inertia.  Hence,  the 
structure  of  the  model  is  generally  composed  of  a  heavy  core  with 
a  light  outer  shell.  The  outer  shell  is  made  of  very  light  poly¬ 
mer  material  and  lead  beads  are  used  as  core  to  regulate  the  posi¬ 
tion  of  the  center  of  gravity  of  the  model.  The  rotational  moment 
of  inertia  obtained  is  (l-2)xl0  3  gm-cm-sec  .  The  geometrical 
dimension  and  the  physical  parameters  of  the  model  must  be  mea¬ 
sured  accurately  before  conducting  the  experiment.  Typical  dimen¬ 
sions  of  the  conical  models  used  in  this  experiment  is  listed  in 
the  following  table  (Table  1): 


TABLE  1.  Typical  geometrical  dimensions  and  physical  parameters 
of  the  conical  models 

[ 


Model 


Geometrical  Dimensions 


Physical  Parameters 


type 

model  code 

semi¬ 

vertex 

angle 

9C  deg 

i 

Rela¬ 
tive 
ctr . 
of 

grav¬ 

ity 

x  /L 
CS 

Rotational 
moment  of 
inertia 

I  gm-cm-sec 

conical 

conical 

ran 

mm 

mm 

Q 

■ 

0.39 

0.41 

1.52xl0-3 

l.llxlO-3 

2 .  Gathering  and  analyzing  the  data 

To  improve  the  precision  in  the  detection  of  the  angular 
position  of  the  model,  the  coordinated  reading  method  is  employed. 
12  points  are  detected  on  the  profile  of  the  model  image  in  the 
HCZ-1  three-dimensional  detector  and  the  probability  error  of  the 
orientation  angle  is  about  ±  0.05°. 
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1 

I 


For  the  plane  free  flight  motion  of  axi-symmetrical  model 
without  rotation,  the  law  of  angular  motion  can  be  written  as 


m 


where 


v- 


Ma 

M 


derivative  of  moments  due  to  aerodynamic  drag, 
if  expressed  in  terms  of  coefficients: 

(Cmq+C")'7--‘3-S 

—  derivative  of  static  pitching  moment  C*»* d*q*S 

—  additional  moment  due  to  slight  asymmetry 


Reference  [2]  has  made  detailed  description  according  to  the 
three-cycles  theory.  Since  under  general  condition. 


(5) 


Hence,  after  linear  assumption  for  angular  motion  of  single  degree 
of  freedom.  Equation  (4)  may  have  solution  of  the  following 
simple  form: 


where 


0  •  K  ’  tu  •  co$(ut  +  9)  +  K, 


(6) 


X 


0:  - 


~cm,  ’  *  •  q  • 

M 


(£»«  +  Cmi )  ~  jjg)  •  Cto]  •  *  »  1  f 


8/  •  V 


(7) 

C8) 


As  a  special  case  of  the  three-cycle  theory,  this  solution 
can  be  represented  by  a  rotating  vector  with  the  pitching  angle  6 
as  the  projection  of  this  vector  on  the  vertical  axis.  The  mag¬ 
nitude  of  the  vector  is  K.  It  rotates  at  an  angular  velocity  of 
oj  and  6  is  the  initial  angular  position.  The  in  Equation  (6) 

is  a  small  adjustment  angle  due  to  the  moment  M  . 

^  a 
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With  Equation  (6)  and  the  data  set  of  angular  motion  of  the 
model  (0,,t.)  ,  ,  -  ,  ...»  taken  from  the  experiment  and  employ- 
ing  the  least  square  method  to  match  them,  the  coefficients  in  the 
equation  can  be  determined. 

The  convergence  criterion  in  the  iterative  matching  process  is 

| SSRJSS Rt  —  1|  <  10“*  (9) 

where 

SSR  —  2  {£  —  t*«  •  **•*'  *  co *(*>,/,  +  *,)  +  X*]}*  ,  . 

<«i  *  ( ID ; 

in  which  the  subscript  o  indicates  the  first  approximation  value 
of  the  parameter,  while  2  and  1  indicate  the  results  obtained  after 
two  successive  iteration.  After  the  iteration  is  completed,  the 
static  and  dynamic  stability  derivatives  of  the  pitching  motion  of 
the  model  can  be  determined  according  to  Equations  (7)  and  (8)  and 
the  u>  and  X  obtained. 


3 •  Results  and  discussion 


Table  2  lists  the  experimental  result  of  static  and  dynamic 
stability  derivatives  for  ©c  *  10°  and  11°  cones  at  =  9.0. 


Table  2.  Experimental  result  of  static  and  dynamic  stability 

derivatives  for  9  =  10°  and  11°  cones 

c 


m 

■ 

m 

m  n  t  1  Jfinff 

motion  parameters 
obtained 

aerodynamic 

derivatives 

matching 

error 

3t  no. 

aynamiL 

pressure 

model 
no . 

.  initial 
angle  of 
attack 

K, 

K 

X 

m 

u>d 

V 

Cmf 

Cm,  +  C»j 

r*j 

a 

1469 

0.182 

10-7 

0.0135 

0.119 

-10.69 

1158 

0.0086 

-1.76 

-3.4 

0.00056 

0.0025 

1470 

0.192 

10-4 

10* 

0.0154 

0.158 

-16.62 

1061 

0.0079 

-1.62 

-6.9 

Q. 00076 

0.0031 

1478 

0.178 

10-3 

12* 

0.0151 

0.169 

-23.00 

1200 

0.0089 

-1.75 

—7,2 

0.00229 

0.0047 

1476 

0.183 

10-6 

16* 

0.0141 

0.260 

-3.64 

977 

0.0072 

-1.53 

-1.2 

0.00181 

0.0042 

1432 

0.186 

11-4 

4* 

0.0117 

0.078 

-9.08 

1249 

0.0091 

—  1 .54 

**  2*2 

0. 00056 

0.0027 

1441 

0.189j 

ll— 9 

12* 

0.0220 

0.119 

-15.36 

1201 

0.0068 

-1.46 

—3.9 

0.00628 

0.0084 

1451 

0.194' 

11-13 

16* 

0.0147 

0.239 

-12.32 

1170 

0.0085 

-1.42 

—3.2 

0.00123 

0.0041 

1471 

0.t95' 

11-6 

10* 

0.0060 

0.111 

-12.79 

1218 

0.0069 

-1.50 

—3,3 

0.00433 

0.0063 

1472 

0.192 

11-7 

20* 

0.0103 

0.267 

-10.79 

1200 

0.0068 

-1.46 

—  2.7 

0.00752 

0.0063 

1474 

0.196 

11-14 

26* 

0.0160 

0.426 

-12.72 

1192 

0.008^ 

-1.55 

-3.3 

0.00199 

0.0047 
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The  table  does  not  only  list  the  motion  parameters  and  the 
aerodynamic  stability  derivative  of  each  experiment,  but  also 
the  probability  error  a  of  matching.  It  reflects  the  deviation 
of  mathematical  model  from  the  actual  angular  motion.  It  is  given 
by  the  following  expression: 


N  is  the  number  of  unknown  coefficients  during  the  matching  pro¬ 
cess.  Here  N  =  5^  n  is  the  number  of  data  points  involved  in 
the  matching. 


As  a  typical  angular  motion  of  the  model,  experiment  1^69  is 
sketched  in  Figure  3-  The  solid  line  in  the  figure  shows  the 
motion  based  on  the  substitution  of  the  gasdynamic  parameters 
obtained  from  the  matching  into  Equation  (6). 


The  static  stability  derivative  C  indicate  that  within  the 
test  range,  the  value  of  Cma  is  independent  of  the  initial  angle 
of  attack.  It  Is  mainly  affected  by  the  location  of  the  center 
of  gravity  of  the  model.  This  shows  a  relatively  good  linearity 
of  the  static  stability  of  the  conical  models.  The  results  are 
plotted  in  Figure  4  and  the  values  are  found  to  be  consistent  with 
Newton’s  theory  which  governs  the  solid  line  in  the  figure  accord¬ 
ing  to  the  following  expression  [33: 

cm%  -  -2.083  (-SjSc.  -  ^  «g  0. .  eo.*,)  (-12) 

The  result  for  the  dynamic  stability  derivative  can  be  seen 
In  Figure  5.  It  is  observed  that  the  pitching  resistance  moment 
coefficient  obtained  is  much  higher  than  the  Newton's  value  and 
that  the  dispersion  is  relatively  large  for  different  experiments. 
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Figure  3.  Typical  detecting  angular  motion  and 
matched  angular  motion 

o  detected  motion  - matched  motion 


M  =9-0,  Re  =1.6x10*  ^  =0.0085,  test  No.  1469,  model  no. 10-7 

C  =-1.76  C  +C  =-3-4 
ma  mq  mot 
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Figure  4.  Experimental  Figure  5.  Experimental  result 

result  of  static  stability  of  dynamic  stability  for 

for  0  =10°,  11°  cones  9  =10°, 11°  cones 

The  dispersion  of  the  experimental  results  is  closely  related 
to  ‘■’"•e  inaccuracy  in  the  determination  of  flew  parameters.  As 
indicated  by  Equation  (8),  the  final  value  of  the  aerodynamic 
resistance  coefficients,  besides  based  on  the  resistance  fa'  ors 
obtained  from  matching,  are  also  determined  by  the  dimension  of 
model  and  the  flow  parameters.  The  determination  of  the  dynamic 
pressure  value  is  especially  difficult  due  to  possible  error  of 
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relative  large  magnitude.  The  general  trend  reflected  by  several 

sets  of  experimental  data  is  much  larger  than  the  Newton's  value. 

The  causes  of  this  phenomena  requires  further  analysis.  From 

the  experimental  result  of  the  effect  of  Reynold's  numbers  and 

the  decrease  in  frequency  on  the  aerodynamic  resistance  given  by 

[It],  we  observe  that  the  dynamic  stability  multiplies  with  the 

increase  in  frequency  at  low  Reynold's  number.  Hence,  with 
4 

ReD  *  3x10  based  on  the  characteristic  length  of  the  model  in 
our  case  and  with  a  low  frequency  f  =  ^  0.01,  pitching  resist¬ 

ance  derivatives  higher  than  the  Newton's  value  are  expected  con¬ 
sequence.  Furthermore,  it  should  be  pointed  out  that  the  models 
used  here  are  not  sealed  at  the  base.  Hence,  there  practically 
exists  a  large  concave  base  which  has  a  considerable  effect  on 
the  aerodynamic  resistance  coefficient.  As  mentioned  before,  the 
dispersion  of  the  resistance  coefficient  is  quite  considerable 
and  the  results  are  essentially  preliminary.  However,  such  dis¬ 
persion  does  not  overshoot  its  order  of  magnitude.  Also,  the  con¬ 
sistency  in  the  characteristic  of  the  dynamic  stability  reflected 
by  different  experiments  demonstrates  the  possible  prospect,  of 
the  measurement  of  dynamic  derivatives  with  the  free  flight  method 
in  the  shock  tunnel. 
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LEAST  SQUARES  FINITE  ELEMENT  ANALYSIS  OF  STEADY  HIGH  SUBSONIC 
PLANE  POTENTIAL  FLOWS 


Jiang  Bo-nan  and  Cbai  Jia-zben 


1.  Method.  The  fundamental  dimensionless  equations  for 
steady,  subsonic  plane  potential  flow  is 


ur  —  V.  —  0,  —  1  +  —  y-1  [A/i  —  C*1  + 


(1) 


where  x  and  y  are  the  orthogonal  coordinates  normalized  by  the 
characteristic  length  of  the  flow  field,  u.  v  and  a  are  respect¬ 
ively  the  velocity  components  and  the  local  sound  speed  normal¬ 
ized  by  the  sound  speed  of  the  free  stream,  M^  is  the  free  stream 
Mach  number  and  y  is  the  specific  heat  ratio. 


The  boundary  conditions  depend  on  the  particular  problem. 

We  will  employ  the  iteration  method  to  find  the  numerical  solution 
of  the  quasi-linear  first  order  partial  differential  equation  Cl). 
At  each  Iteration  cycle,  the  least  square  finite  element  method  is 
used  to  solve  the  linearized  equation  set.  Both  the  variable  stiff 
ness  and  the  constant  stiffness  models  are  tried  out.  Their  diff¬ 
erence  is  reflected  in  the  linearized  equations,  that  is,  the  coeff 
icients  are  variables  with  the  former  model  but  are  constants  with 
the  latter. 


1)  Variable  stiffness  method.  Let 


Then  Equation  (1)  can  be  rewritten  as 

iu,  +  Bp,  +  C(n,  +  *7  *#  ■»  0  (3) 
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Coefficients  A,  B,  C  are  regarded  as  known  quantities  and 
are  determined  by  the  velocity  field  solved  in  the  preceding  iter¬ 
ation  cycle.  The  new  values  of  u  and  v  are  determined  with  the 
following  method.  Firstly,  the  flow  field  are  divided  into  finite 
elements.  In  order  to  simulate  boundaries  of  complicated  pro¬ 
file,  quadrilateral  elements  with  eight  nodal  points  are  employed 
[1].  The  two  velocity  components  at  each  node  are  taken  as  the 
main  unknowns.  The  coordinates  and  the  unknowns  of  any  point  in 
the  element  can  be  expressed  as 


(4) 

r.K*  **•<} 

(5) 

where  N^CSjn)  is  the  shape  function,  x.^,  y.^  are  the  coordinates 
of  the  i  th  node  of  the  element,  u^,  are  the  unknown  velocity 
components  of  the  i  th  node  of  the  element . 

Equation  (5)  is  practically  the  expression  of  approximate  solu¬ 
tion  of  equation  (1).  Generally,  Equation  (1)  are  not  satisfied 
exactly  no  matter  how  u.^  and  v.^  are  chosen.  That  is  for  any  point 
in  the  element,  there  exist  residues: 

*,  -  2  IGWfc.  +  CN,.>,  +  +  CN,.M  ... 

i-i  t  (o; 

#•1 

Employing  the  least  square  method  [2],  u^^  and  are  chosen  such 
that  the  residues  are  minimized,  that  is  to  construct  a  second  order 
functional  (assuming  that  there  is  only  one  element) 

/-JJ  (*}+•*})**)-  (7) 

Received  October  IT,  1978. 
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In  the  above  expression,  the  integration  is  conducted  in 
the  region  Sg  and  a  is  a  positive  constant.  Generally,  a  =  1  is 
taken.  This  indicates  that  the  continuity  equation  and  the  irro- 
tational  equation  are  equally  important.  Taking  very  small  value 
of  I  with  respect  to  ,  v^,  the  linear  algebraic  equations  with 
respect  to  v.^  are  obtained: 

(P 

[K„l  <#!*>■ 

I  Kh]  [KB1 

•  •  • 

•  #  * 

•  •  •  * 

I**)  (!?■)•' •M«]J 

Key:  (1)  cemetery 

In  which  the  sub-matrix  of  the  stiffness  matrix  of  the  element  is 

■AtNt„Ni..+  AC(.N„M„  '  !  A  +  ACN,..N)., 

+  +  (C>  +  a)NllfNhr  [  +  9Cff,.,N,.f  +(C*-a)Nt.,N,.. 

jj  . . . . 

*  ABNt.JI,.'  +  ACN,.Jt,.'  |  +  BC(.Nt„N,., 

+  B CN^N,.,  +(C,-«)W,^„  :  +  +  C C*  + 

The  integral  In  the  above  equation  Is  calculated  numerically 
by  2x2  Gaussian  integration  formula.  The  sub-vector  (F^}  of  the 
loading  vector  of  the  element  is  the  zero  vector. 

For  aggregates  of  many  elements  the  overall  rigidity  matrix 
and  the  overall  load  vector  can  be  obtained  on  the  basis  of 
normal  aggregate  rules.  After  considering  boundary  condition 
corrections,  the  entire  set  of  equations  can  be  obtained.  It 
should  be  noted  that  the  entire  load  vector  which  has  undergone 
boundary  condition  correction  is  no  longer  a  zero  vector. 

2)  Constant  stiffness  method. 

> 

Equation  (1)  is  rewritten  as 

#«  +  •>»/,«,-'  r,  m  0  (1°) 

where 

/  —  f  I  +  J-~  1  (W»  —  (**  +  +  *r(*,  +  r.)  +  •’*»'*] 

12  ‘  (21) 
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Here  f  is  known,  and  is  determined  by  the  velocity  field  given 
in  the  preceding  iteration.  Taking  A  *  B  *  1,  C  ®  0  in  Equation 
(9),  the  expression  for  the  stiffness  sub-matrix  of  the  constant 
stiffness  method  is  obtained.  The  sub-vector  of  the  loading  vec¬ 
tor  of  the  element  is 


(12) 


2 .  Some  problems  in  the  calculation. 

Since  the  velocity  component  normal  to  the  solid  boundary  is 
zero,  the  inclined  boundary  condition  is  utilized,  that  is,  at  the 
nodal  points  on  the  solid  boundary,  the  tangential  and  the  normal 
velocity  components  are  taken  as  unknowns. 

The  symmetry  of  the  stiffness  matrix  is  favorable  to  the  cal¬ 
culation.  The  following  methods  are  used  to  solve  the  linear  alge¬ 
braic  equation:  the  triangular  decomposing  method  of  variable  band 
width  one-dimensionally  stored  total  stiffness  matrix  and  the  elim¬ 
ination  method. 

Set  f  =  0  in  Equations  (10),  and  the  flow  field  obtained  by 
solving  this  incompressible  problem  is  taken  as  the  initial  field. 
At  all  the  nodal  points  of  the  field,  the  problem  is  considered 
convergent  if  the  relative  variation  of  two  successive  values  of 
the  dimensionless  density  p  *■  [l  +  —  W  +  r*))J 

is  less  than  a  given  small  quantity  e. 

3 .  Numerical  examples 

Several  typical  examples  are  calculated  with  ALGOL  60  language 
on  the  TQ-16  machine. 


it  i 
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1)  Plow  past  a  circular  cylinder.  The  flow  condition  at 
seven  times  the  radius  of  the  cylinder  is  considered  undisturbed 
from  the  free  stream.  One  quarter  of  the  circular  cylinder  is 
divided  into  2x2  or  4x5  elements,  e  is  taken  as  0.0001  and  0.0009 
respectively.  The  calculated  results  are  shown  in  Table  1  and 
Figure  1.  It  can  be  observed  that  even  when  reaches  the  crit¬ 
ical  Mach  number  0.42,  the  variable  stiffness  method  is  still  con¬ 
verging  rapidly. 


Table  1.  Number  of  iterations  for  flow  past  circular  cylinder 


X  No-  0f  \  M 

CalculationX.  inter-  X  “ 


scheme 


20 

elements 


\ations  \ 


ti 

d 


variable  stiff- 
ess  method 


constant  stiff- 
neww  method 
variable  stiff¬ 
ness  method 


0.20 

o 

cn 

o 

O 

-tT 

o 

2 

3 

8 

2 

3 

■1 

3 

6 

27 

2 

2 

3 

8  32  no  convergence 


27  I  no  convergence 


h  a 
JS' 


,V«.-0.49  . 

•  V 


Figure  1.  Mach  number  on  the  surface  of  the  circular 
cylinder 

o  finite  element  (variable  stiffness  method) 

- result  from  reference  [33 


2)  Plow  past  a  symmetric  airfoil.  The  symmetrical  NACA 
0012  airfoil  with  zero  angle  of  attack  is  calculated.  The  upper 
half  region  of  the  airfoil  is  divided  into  5x12  elements  and 
e  =  0.0001  is  taken.  Convergence  is  achieved  with  the  variable 
stiffness  method  after  six  iterations  (Figure  2).  There  is  no 
convergence  if  the  constant  stiffness  method  is  used. 

4 .  Discussion 

With  the  above  presentation,  it  can  be  concluded  that  the 
least  square  finite  element  method  is  feasible  for  steady  sub- 
onic  potential  flows.  The  variable  stiffness  method  is  especially 
suitable  for  high  subsonic  flows. 


Figure  2.  Mach  number  on  the  airfoil  surface 
o  finite  element  (variable  stiffness  method) 

-  result  from  reference  [4] 

NACA  0012  M^  =  0.72  0°  angle  of  attack 

Nowadays  many  finite  element  methods  [5]  have  been  proposed 
for  solving  problems  of  compressible  potential  flows.  The  stream 
function  ip  or  the  potential  function  <p  is  considered  as  the  main 
unknown  quantities  in  most  methods.  There  are  several  advantages 
of  the  present  method.  Firstly,  the  velocities  at  the  nodes  are 
obtained  directly,  avoiding  any  numerical  differentiation.  Hence, 
relatively  coarse  mesh  can  be  employed.  Secondly,  the  difficulty 


of  determining  the  density  at  high  subsonic  Mach  numbers  when 
using  is  also  avoided.  Thirdly,  the  boundary  conditions  are 
simpler,  in  contrast  with  the  cases  when  b  is  employed.  Fourthly, 
the  continuity  requirement  of  the  shape  function  is  reduced  from 
C’  to  C°  as  compared  with  the  least  square  method  for  4>  or  <p . 

Finally,  we  thank  Professor  Lo  Shi-jun  for  his  valuable  advice 
and  Mr.  Sun  Huei-min  for  revising  the  computer  program. 
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THE  AERODYNAMICAL  ANALYSIS  OF  BODY  IN  HYPERSONIC  SOURCE  FLOW 
FIELD 


Ling  Guo-can 

(Institute  of  Mechanics,  Academia  Sinica) 

The  hypersonic  source  flow  field  is  a  typical  non-uniform 
free  stream  flow  field.  Many  high  speed,  high  enthalpy  experi¬ 
mental  systems,  such  as  shock  tunnels  equipped  with  conical 
nozzles,  the  gun  wind  tunnels,  the  electric  arc  wind  tunnels, 
etc.,  are  characterized  by  hypersonic  spherical  expansion.  Cer¬ 
tain  external  flows  around  bodies  also  display  the  same  character¬ 
istics.  Simple  but  accurate  analytical  methods  are  required  to 
determine  the  aerodynamic  force  on  the  body  and  the  effects  on  the 
pressure  distribution  and  various  aerodynamic  characteristics  for 
such  non-uniform  free  stream  condition.  Some  efforts  have  been 
made  [1-7]  but  no  satisfactory  analytical  results  have  been  pre¬ 
sented  so  far.  Based  on  the  characteristics  of  the  hypersonic 
source  flow  field  and  the  simple  form  of  the  potential  function, 
the  present  study  obtained  various  aerodynamic  coefficients  by 
means  of  Newton's  theory.  Under  the  condition  of  small  angle  of 
attack,  the  analytical  results  are  satisfactory.  Due  to  practical 
requirements,  analysis  on  the  effects  of  the  non-uniformity  of  the 
free  stream  is  emphasized.  Agreement  between  the  calculated  results 
and  the  experimental  data  is  satisfactory. 

1.  Hypersonic  source  flow  field 

p 

At  hypersonic  conditons,  M>>1,  and  when  M  (y-1)/2>>1,  the 
following  relations  among  physical  properties  in  steady,  ideal 
point  source  flow  are  obtained  by  integrating  the  basic  equations: 

V  *  constant,  P  -  P  -  (hCr/r,r\  U  (1) 

when  r  is  the  distance  of  any  point  in  space  from  the  point  source. 


The  subscript  denotes  certain  reference  values  in  the  flow  field. 
Note  that  the  velocity  of  the  above  flow  field  is  almost  constant . 
Hence,  the  corresponding  potential  function  varies  with  r  only, 
that  is  • +  rV  +  C  .  This  potential  function  provides  some  con¬ 

venience  for  the  analysis  following.  If  the  point  source  is 
located  at  the  coordinates  (f»  0,  C),  on  the  symmetric  surface  xoz  of 
the  Cartesian  coordinates  Oxyz,  then  the  velocity  potential  and 
the  velocity  distribution  should  be: 


9  -  vVc*  +  sy  +  /  +  c*  +  c  y  +  c 

V,  -  (x  +  L.co*«)KL-*,  V,-yVL'\  V, .  -  («  +  L,wa>VLr' 
L  *■>  t/Cx  +  L0co»a)’  +  y*  +  C*  +  Lgiina)' 


(2) 


in  which  Lq  is  the  distance  from  the  point  source  to  the  origin, 
and  a  is  the  angle  between  Lq  and  the  x-axis . 


2.  Pressure  distribution 


Let  the  geometry  of  the  body  of  revolution  in  the  flow  field 
described  above  be  ,  and  the  unit  vector  normal  to  the 

body  surface  be  n*  —  «m4l  —  co»aco»0j  —  eottunfik,  a  — 

The  elevation  angle  of  the  line  Joining  the  point  source  and  the 

vortex  of  the  body  from  the  Ox-axis  is  defined  as  the  angle  of 

attack  a  of  the  non-uniform  free-stream.  The  Newton's  hypersonic 

impulse  theory  gives  the  pressure  coefficient,  C„  —  p*(V»  • 

of  any  point  on  the  upwind  surface  of  the  body.  The  subscript  b 

denotes  the  local  freestream  parameters  at  any  point  on  the  body 

surface  given  by  Equations  (1)  and  (2),  and  x,y,z  satisfy  the  equa- 

2 

tion  of  the  body  profile.  qQ  »  pQV  /2  is  the  local  free  stream 
dynamic  pressure  at  the  vertex  of  the  body.  If  the  distance  Lq 
from  the  point  source  to  the  body  vertex  is  taken  as  the  character¬ 
istic  length  of  the  flow  and  the  reference  value  for  normalization, 
then  the  pressure  coefficient  at  any  point  on  the  upwind  surface  of 
the  body  is 
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Apparently,  when  the  point  source  is  infinitely  far  away  from  the 
body,  the  above  expression  reduces  down  to  the  well-known  Newton’s 
formula  of  the  uniform  flow. 


For  conical  bodies  f  — *tg 8,  and  the  pressure  distribution 
on  the  conical  surface  is 

Cfi  ■*  C,UL~*,  CpU  —  21  cot  a  tin  3  —  tin  a  co«  6  lin  0]’  1 
L  1 1  4  P  +  2 £ coso  +  +  2*tg  Jiinatia/fji  i 

The  subscript  U  denotes  the  uniform  free  stream  values.  Hence, 
there  exists  a  simple  transformation  relation  between  the  pressure 
coefficient  on  the  conical  surface  in  source  flow  field  and  the 
pressure  coefficient  in  the  corresponding  uniform  flow  field.  The 
former  can  be  expressed  as  a  function  of  the  uniform  flow  values 
and  L,  the  dimensionless  distance  from  the  point  source  to  the 
conical  surface.  The  pressure  coefficient  decreases  along  the 
conical  surface,  hence  losing  the  characteristics  of  conical  flows. 
This  effect  is  greater  when  the  point  source  is  closer  to  the  body. 

For  a  spherical  surface  of  radius  p  —  R‘— (R— i)1,  .  At 

zero  angle  of  attack,  the  pressure  coefficient  on  the  upwind  sur¬ 
face  is 

C  -  Q+  dag]1 

"  “  (1  +  2RC1  -  naO')  -  *»*)?  (5) 

where  8  is  the  circumferential  angle  as  shown  in  Figure  2 .  At  the 
stagnation  point,  9  *  it/2.  The  region  of  action  is 

Figures  1  and  2  present  the  calculated  relative  values 
(CrB  —  Cf,}/C„  and  of  the  pressure  distribution  on  the 

conical  surface  and  the  spherical  surface  respectively. 

This  paper  was  received  on  December  18,  1978. 

1)  This  paper  was  presented  at  the  Chinese  Institute  of  Mechanics'  First  Convention 
of  Shock  Tunnel,  November  4-12,  1978. 
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Figure  1.  Pressure  distribution  on  the  conical  surface  (a=0°) 
- calculated  values  according  to  Equation  (4) 


A  m-io.ii  a  — 14*  gun  wind  tunnel  experiment  [5] 
o  «-9.s  #-is*  shock  tunnel  experiment  (large  cone)  [3] 

•  *-9.s  a  —  is*  shock  tunnel  experiment  (small  cone)  [3] 

«  w-a.7  a  -  to*  shock  tunnel  experiment 

■  *-7.3  #-i3*  shock  tunnel  experiment  (large  cone)  [3] 

o  *-7.3  a  —  is*  shock  tunnel  experiment  (small  cone)  [3] 


Figure  2.  Pressure  distribution  on  the  spherical  surface  (a=0°) 
- calculated  values  according  to  Equation  (5) 

•  X - o.oo9  *-i5.9  sphere-cylinder  experiment  [5] 

♦  *-0.0337  v  -  i.e  sphere-cone  experiment 

a  X- 0.0303  *-io.3i  sphere-cone  experiment  [9] 


3 .  Aerodynamic  coefficients  of  typical  bodies 


By  applying  the  coefficient  formula  (3)  obtained,  the  nor¬ 
mal  force  coefficient,  axial  force  coefficient,  pitching  moment 
coefficient  and  the  pressure  center  coefficient  of  the  body  of 
revolution  at  any  angle  of  attack  can  be  calculated  directly: 


c„, - ^  F'  (*'  c,j  un  fidfidi,  C1t  -  4t  (*’  \*'  C'ffdfidx 

cm, - -  ('•  ('*  C,Ki  +  Tnunpdpdx,  ic  -  Cu,lCH. 

xry  '* 


where  ;c>  —  xCf/l > *  and  rg  are  the  body  length  and  the  base  radius 
respectively.  The  upper  limits  of  the  integrations  *»>  are 

the  dividing  line  of  the  upwind  and  downwind  surfaces  and  are 
determined  by  the  zero  pressure  condition  on  the  body  surface. 

For  conical  bodies,  a  <  -  «/2.  f.  -  «in-'Ogtf/»g«),  *,  - 

For  spherical  bodies ,  -  *,-{[« «««-<*  +  -  «)']-*< *ino>“‘h 

For  small  angle  of  attack  when  I/j?  <  (1  +  R)**1  fl,  —  m/2. 

For  small  angle  of  attack,  the  higher  order  terms  of  a  can  be 
neglected.  Also  when  the  point  source  is  considerably  far  from 
the  body,  the  distance  of  any  point  on  the  body  from  the  point 
source  is  almost  equal  to  its  projection  on  the  x-axis.  Simple 
analytical  expressions  for  aerodynamic  coefficients  are  obtained 
through  integration.  The  position  of  the  point  source  still  has 
a  sufficiently  large  range  in  these  expressions.  For  cones  at 
a<9  ,  there  are  certain  correlations  between  the  aerodynamic 

coefficients  in  the  point  source  flow  field  and  the  corresponding 
coefficients  in  the  uniform  flow.  The  normal  force  coefficients, 
axial  force  coefficients,  static  derivative  pressure  center 

coefficient,  pitching  moment  coefficient  in  the  point  source  flow 
field  can  all  be  expressed  by  the  product  of  the  corresponding 
value  in  the  uniform  flow  field  and  the  correlation  function 
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(7) 


Cn,  *“  Cvt/<PiCO»  Cj,  Cni^iCOr  Gm—  m  Cfi,a<p£l')  1 

*C„  -  *Ctu<PlO'),  Cms  -  CmU<Pi(.OfiO)  J 

CNU  —  20C0*1 Ctv  ™  C«**  +  2:^  ff)eo»J  C*m  “  2co*J  A 

*«»  “  2/3  ,ec,*»  c-»  ■“  4/3  « 

«.,(/)- C3  +  o/3< i +  /y,  ^/)^p  +  i/3/r  ^ 

The  correlation  functions  <pXO  depends  on  the  dimensionless 

characteristic  quantity  jt  only,  which  denotes  the  relative  posi¬ 
tion  of  the  body  in  the  point  source  flow  field.  Figure  3  shows 
the  variation  curve  of  1  —  versus  !  .  The  relative  variation 
of  the  aerodynamic  coefficients  due  to  the  non-uniformity  of  the 
free  stream  under  the  condition  of  different  point  source  posi¬ 
tions  can  be  determined  according  to  this  curve. 


SYMBOL  CONTENT 


M  Ref 


• 

pressure  center  experiment 

0 

7 

small 

16 

♦ 

pressure  center  experiment 

0 

6 

small 

20 

5 

■ 

pressure  center  experiment 

0 

10 

7,5,10 

8.2 

o 

pressure  center  modification  calculation 

0  small 

small 

5 

A 

pressure  center  modification  calculation 

16 

2-15 

7 

V 

normal  force  coefficient  experiment 

D 

7 

5 

16 

5 

0 

axial  farce  coefficient  modification  calculation 

b  .006 

16 

2-15 

7 

V 

normal  force  coefficient  modification  calculation 

0.006 

16 

2-15 

7 

For  the  spherical  section  and  for  J. <(1 +  £)'"’  ,  the  following 

expressions  can  be  obtained  after  some  complicated  integrations: 


Cut  * 


CJt  » 


-uM  *  «• + + + *»•' 

R  R  1 

+ 1  +  2*0  +  *y.  +  fnri*oVl)0‘f 

mnrmr^u.  f5' + «• + « + *«•  ♦  *»  +  *xj 


+ _ f  1  ....} 

1  +  2*0  +  *)V 


ChJU 


(10) 


where 


*,  -  2.5  +  14*  +  28*'  +  24**  +  8*«f*  *,  -  2*Q  +  3*  +  2**) 

R>  -  -C3  +  12*  +  lfi*,+8*,)t  *,  -  -C3  +  18*+40*,+40*,+  l6*‘) 
*,  -  0.5  +  4*  +  12*>  +  16**  +  8** 

S,  -  l+7*+2Q*,+30*’+24*«+8*»,  $,r4*+20*,+32*,+20*«+4*5 
S,  -  -(2**  +  «**  +  «*♦  +  2*’),  S,  -  -<3+iS*+28*,+24*,+8*0 
S,  —  -O  +  7*  +  20  *‘  +  30*’  +  24*4  4-  8*0  i: ;V.t' "  ' 

*  -  */U,  1.-I./R  ' 


66 


/.  is  the  thickness  of  the  spherical  body  in  the  x-direction. 
When  i?-»o  j  the  above  limits  become  the  Newton's  formulae  in  the 
uniform  flow.  Figure  4  shows  the  variation  curve  of  the  deriva¬ 
tives  of  normal  force  coefficient  and  axial  force  coefficients  of 
the  spherical  section  vs.  the  point  source  position. 


0.8OHU0 


Figure  4 .  and  CTs  of  the  spherical  section 

For  a  slender  cone-sphere  combination  of  bluntness  n  at  asS, 
and  for  sin  6^6,  cos  <5^1,  the  following  can  be  obtained: 

Cmmu  “■  ij*  4-  C tit,,  Cjuau  “  Curf  +  Cn,,  Cum m„if  +  CK.i,  1  -  . 


*c»— [C«‘ \-l(\  +  C"*']/Cw“* 


Figures  1-3  and  Figures  5  and  6  show  the  comparison  between 
calculation  results  and  the  experimental  data.  The  force  mea¬ 
surement  accuracy  of  current  impulsive  wind  tunnel  (about  5-10?) 
is  one  order  of  magnitude  lower  than  that  of  the  ordinary  wind 
tunnel.  Some  of  them  reach  around  17?.  Since  there  is  a  lack  of 
experimental  data  for  single  spherical  section.  Figure  2  adopts 
the  results  from  the  sphere-cone  and  the  sphere-cylinder.  Hence, 
the  data  close  to  the  stagnation  point  at  the  head  should  be  taken. 
The  accurate  solution  of  the  pressure  distribution  on  the  cone  of 
the  uniform  flow  is  taken  from  [10].  The  pressure  distribution  and 
the  relative  variation  of  aerodynamic  coefficients  due  to  hyper¬ 
sonic  point-  source  free  stream  calculated  here  are  found  to  match 
with  the  experimental  results.  The  calculated  aerodynamic  coeffi¬ 
cients  of  the  sphere-cone  combination  also  agree  quite  well  with 
the  experimental  results.  They  are  also  consistent  with  the  modi¬ 
fied  values  given  by  [5,7].  Compared  with  the  uniform  free  stream 
condition,  the  source  flow  causes  the  pressure  center  of  the  sphere-cone  to  move 
forward  and  the  normal  and  axial  force. goefficients  to  decrease.  For  slender 


Figure  5.  Normal  force  coefficient  and  pitching  moment  coeffi¬ 
cient  of  sphere-cone 
4  calculation  by  Equations  (12),  (13) 

•  gun-wind  tunnel  experiment  [7]  <$  =12.5°  n  =0-0.5  a=2-l 5° 

a  gun  wind  tunnel  experiment  [7]  5  =16°  n  =0-0.5  a=0-l 6° 

*  gun  wind  tunnel  experiment  [7]  5  =20°  n  =0-0.5  a=0-17° 

-o  uniform  flow  values  modified  by  the  present  method 

o  uniform  flow  values  modified  by  [7] 
a  uniform  flew  values  modified  by  [7] 

0  uniform  flew  values  modified  by  [7] 

- experimental  values  of  uniform  flow  [12] 


Figure  6.  Pressure  center  of  sphere-cone 
k  shock  tunnel  experimental  values 

A  pressure  center  of  the  uniform  flow  modified  by  the  present 
method 

n-0.15  ct=10°  M-8.2,  13.5  n-0.15  a=10°  M-8.2,  13-5 

- calculated  value  according  to  the  uniform  flow  theory 

n*0.15  oc*10° 
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cones  at  small  angle  of  attack,  the  pressure  center  moves  forward  by  1— ?//). 
When  J>0.O3 ,  that  variation  is  more  than  1%.  The  relative 

variation  of  both  the  normal  and  axial  force  coefficients 
is  1— 9»<0  .  When  />0.04  a  the  variation  is  more  than  10.%. 

For  the  hypersonic  uniform  free  stream  problem,  the  Newton's  theory 
is  effective  in  determining  the  pressure  distribution  and  aero¬ 
dynamic  coefficients,  provided  that  the  shock  stays  close  to  the 
body  surface  [11,12].  For  hypersonic  source  flow  field,  it  also 
displays  a  satisfactory  accuracy  in  analyzing  the  aerodynamic 
effects  due  to  the  non-uniformity  of  free  stream. 
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